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ABSTRACT 


Porous media are interconnected networks of void 
space, of multitude shapes and sizes. Attempts in developing 
an analytical expression for the fluid flow behaviour through 
them have not been satisfactory to-date. The classical model of 
the bundle of tubes is too simplified a model to be realistic 
and useful. In the present work the flow behaviour of porous 
media has been modelled by a network of interconnected tubes 
of different sizes. An attempt has been made to ext’end the 
concept of relating pore-size distribution (both volumetric 
and number) to the capillary character of the network. The 
validity of ’ink-bottle* effect has been tested by the 
comparison of capillary pressure data with the actual tube size 
distribution of the network. 

The extent of inter connection between the pores in a 
network is found to be a very important factor influencing the 
shape of the P c curves. Therefore, any model of porous media 
based on bundle of tubes concept is likely to give erroneous 
results. Patt’s model, however, may be used with modified value 
of the exponent (between -2.5 to -3.5) to give better result 
for the number distribution. Similarly while correcting the 
measured pore-size distribution (volumetric) by Meyer’s method, 
it has been found that the relation improves the 

interpretation. These values for exponents have been obtained 



(viii) 


for networks withw six tubes (on the average) meeting 

* # 

at a node point. The modified values of exponents perhaps 
correct for interconnection between the pores. 




CHAPTER 1 


INTRODUCTION \ 

mm ,~m- irr run itm-tt. V 

. Porous media, whether in the form of natural rocks or 
particulate aggregates such as powder-packs, have quite a few 
interesting properties compared to the properties of their 
bulk solid phase. One such property, the flow of fluids through 
them, has been engaging the attention of scientists in many 
disciplines. The petroleum engineers are interested in knowing 
the quantity and the maximum possible rate of fluid withdrawal 
from reservoir rocks. The so il-3cientists concern themselves 
with the distribution and movement of soil-moisture* Chemical 
engineers, likewise, dealing with the problem of drying, are 
faced with the question of determining the rate at which the 
moisture from the pores can move up to the surface. Similarly 
while analyzing the performance of a packed bed, it is pertinent 
to evaluate the effectiveness of mass-transfer as well as liquid 
hold-up. Internal void structure is important in determining 
the effectiveness of batteries as well as the moisture movement 
capacity of the concrete and bricks. 

The total internal surface area, porosity or permeability 
do not provide adequate information about the complex pore 
structure of porous media. Quite often, it is necessary not 
only to know these cumulative properties but also the distri- 
bution of pore sizes. The shapes of the relative permeability 
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and. capillary pressure curves are markedly affected by pore- 
size distribution. Of the two, relative permeability does 
not seem to be suitable for determining the pore-size distri- 
bution because of the averaging effect associated with the 
flow phenomena in the whole spectrum of pore-sizes. Capillary 
pressure, on the other hand, is a unique function of pore-size, 
interfacial tension and wettability. Keeping the latter two 
quantities constant it is possible to estimate the sizes and 
proportions of the contributing pores from the successive higher 
values of capillary pressure. 

Capillary phenomena comes into play whenever a fluid/ 
fluid interface is in contact with a solid surface. Due to 
unequal interfacial tensions, fluid adjusts itself to attain 
thermo dynamic equilibrium, thereby min ini zi ng surface free 
energy. The well known Laplaco equation which relates the 
various quantities affecting capillary pressure is 



where = Capillary pressure, 

V 

- Interfacial tension 

r-j&r 2 = Two principal radii of curvature of the 
surface at any point. 

The implicit assumption that the sum of reciprocals 
of the two principal radii of curvature (mean curvature) is 
the same at every point on the equilibrium surface, is a 
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mathematical consequence of the fact that the surface area 

* 

tends to a minimum (i.e, spherical or nearly so in most cases). 
All methods of pore-size determination, that make use of 
capillary pressure are based on the hypothesis that the mean 
curvature of an equilibrium liquid surface can somehow be 
related to the pore-size. This relationship is difficult to 
establish for any shapes other than that with the simplest 
geometry. For straight cylindrical capillaries, Laplace's 


equation can be written as 

n 2 *7 cos & 


( 1 . 2 ) 


c r 

where © is the angle of contact of the interface through the 
wetting phase and r is the radius of the capillary. This 
equation can correctly predict the pore-sizes and their distri- 
bution if bundle of parallel tubes is a valid model for porous 
media. For real system, this is an obvious over simplification 
and the equation will relate the capillary pressure to the 
radius of the largest opening leading to the pore rather than 
the radius of the pore itself. But in the absence of a better 
known method this pore entry radius is taken as a measure of 
the actual pore-size. Implied in this is the assumption that 
the actual pore-size is some unique function of the pore entry 
size for pores of all sizes. This however is difficult to 
justify. 

In capillary pressure experiment with porous media, 
the sample is gradually desaturated of the wetting phase by 
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injecting some non-wetting fluid (making known contact angle 


with the solid surface) at successive higher pressure. At 
equilibrium, the capillary pressure opposing the entry of 
non-wetting phase, P j is balanced by the externally applied 
pressure. This P 0 when put in Eq.(l.2) gives the radius of 
the smallest pore, r, invaded by the non-wetting phase. Knowing 
the saturation of the non-wetting phase in the sample at the 
stage, we arrive at the value of the fractional pore volume 
contained in the pores with radius greater than or equal to r. 

By successively increasing P , non-wetting phase is forced into 
smaller and smaller pores and the fractional pore volume 
corresponding to different values of r can be obtained. This 
method has been universally used for the calculation of volumetric 
pore-size distribution in porous media. 

For fluid flow calculations, the pore number distribution 

is the desired quantity. To arrive at the number distribution 

from volume distribution, it is essential that we make some 

assumption about the pore-geometry and its dimensions. Por 

parallel tube model of porous media, the quantity in question 

is the relation between the tube radius and length. Based, on 

the assumption that 1 = Ch*- , where C is proportionality constant 

1 

and r is radius, Fatt derives the number distribution function 
as follows? s ? 

r -2+^c. . 

J p c &s 

*->-j ^ 



i ( i?2 ) 


( 1 . 3 ) 
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where f(r 2 -r 1 ) = frequency of occurrence of pores in the 

interval r^-r^ 

g = saturation of the non-wetting phase 

There is considerable disagreement between different 

workers in the field about the value that should assume. 

1 2 
Whereas Patt proposes a value of -1. Meyer implicitly takes 

this as 1 . Dallavalle^ seems to have some experimental results 

to back his suggestion of the value between 0 and 1. In the 

absence of a conclusive evidence, this factor still remains 

uncertain. 

The method just discussed is essentially empirical in 
nature, for quantitative description of flow phenomena in 
porous media, a functional relation between the pore-size 
distribution and capillary pressure curve is needed. It will 
be our endeavour in the present work to exanine the existing 
models and suggest one which can reasonably interpret the 
experimenta.1 vs saturation curves in terms of pore-size 
distribution. 


***** 



CHAPTER 2 


LITERATURE REVIEW MD ANALYSIS 

Natural porous systems can be viewed as inter-connected 

network of voids of different shapes and sizes separated by 

interconnecting links or necks. Often the behaviour of such 

a system is approximated by a bundle of capillaries^. Dullien 
5 

and Batra have collected a comprehensive list of references 
of the attempts made in this direction in their review article. 

It is apparent that very few real porous media will contain 
pores of the straight cylindrical shapes, and therefore the 
pore radii calculated from the capillary pressure have to be 
considered as equivalent pore radii. This is defined as the 
radius of a straight cylindrical capillary that would give rise 
to the same capillary pressure as the measured value. Since a 
variety of different capillary shapes may give rise to the same 
capillary pressure, the pore dimensions calculated from capillary 
pressure measurements are at best semi-quantitative. The equation 
relating properties of porous media to the radius distribution 
of the equivalent bundle of tubes have been given by a number 
of authors . However, apart from the mathematical simplicity 
associated, the most glaring weakness of the model is the 
prediction of highly anisotropic properties. Porous systems 
on the other hand are known to be isotropic in nature. Inspite 
of its success in correlating certain properties viz. permeability 
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formation factor etc., the model has failed to approximate 
the capillary pressure vs saturation and its hysteresis. 

for more complex systems, such as random packing of 
non-spherical particles, Carman^ suggests that twice the 
hydraulic radius is the suitable value of r in the following 
expression for P : 


2 Y cos Q 
r 


( 2 . 1 ) 


where the hydraulic radius is defined as the ratio of porosity 

f to the surface area A in unit bulk volume. A more accurate 

account of the curvature for different geometry of the capillary 

”1 o 

is given by Gregg and Sing . On substitution Eg. (2.1) is 
reduced to 


P c = -f S QS-g-i : (2.2) 

Substituting the expression for the specific surface area 

1 1 

in terms of permeability K and porosity, from the Kozeny 
equation, Eg. (2,2) becomes 


_ ( 1 ) 7 * 


1-i 


( 2 . 3 ) 


where the threshold pressure P^ has been substituted for P c . 
k is an empirical textural constant of the medium and accounts 
for the deviation of a porous system from the model. 

This was obviously a gross simplification, and Carman 
himself points out that, for porous bodies with widely varying 
pore-size distribution, calculations based on Kozeny-Carman 
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model can be quite wrong. 


Sphere Pack Model s 


The earliest attempt in approximating the flow behaviour 
of porous media was their simulation by packs of uniform spheres. 
But the complexity of the pore-geometry in sphere-packs prevented 
the derivation of an accurate and meaningful description of 

flow behaviour through them. Recognizing the difficulty, 

1 1 

Koseny inductively developed an equation relating permeability 
to porosity and internal surface area of a sphere-pack. 


where 


K = 


C f 3 
TET 


f = porosity 


(2.4) 


T = tortuosity 
S = specific surface area 

12 

Kozeny equation, as modified by Carman and given as 


K 


• l3 

k S (1-fT 


(2.5) 


did not prove to be of much value, aside from its use in 
estimation of the surface area of powders. 

A comprehensive review of the use of sphere-packs to 

model the flow behaviour of packed beds, has been given by 

1 5 

Haughey and Beveridge . In their analysis, the emphasis has 
been in arriving at the mean fractional voidage for both 
regular and random packing, on the basis of the position coordi- 
nates and the coordination numbers, rather than the pore _ 



9 


geometry. Generally speaking, particle shape and. size 
distribution are the two factors most likely to affect the 
packing structure and its properties. Some work has been 
carried out on the densest packing of non-spherical particles. 
But the multiplicity of packing types with the same coordination 
number and vice versa, accompanied by the size distribution of 
particles further complicate the pore structure. Hence the 
task of analytical description is rendered almost hopeless. 

The problem of capillary hysteresis in uniform sphere 

14 - 15 

packs was studied by ICruyer , Erevel and Kressley , Meyer and 

1 6 1 7 18 

Stowe ’ and more recently, by Melrose . Their analysis 

points out that the curvature of the invading interface of 

non-wetting phase is generally larger resulting in higher value 

of capillary pressure compared to the retracting interface 

associated with imbibition. This, according to them, may be 

the reason for the porous media exhibiting capillary hysteresis. 

IQ 

Naar and Wygal " view any porous media as a random 
mixture of spherical particles and propose a model based on the 
principle of averaging the properties of the basic sets to 
arrive at the mixture-properties. The problem essentially is 
finding of appropriate weighting factors to average the 
basic set properties, determined experimentally. Expressions 
for porosity, permeability and capillary pressure behaviour 
of the porous media have been developed. But the experimental 
support to this model is too meagre to be accepted. 
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Network Model s 

The description of single and multiphase fluid flow 
through capillary pores in porous media has been unsatisfactory 
because the pore-geometry is too complex to be described 
adequately by analytical expression. Markin used a model 
consisting of a system of randomly arranged intersecting voids 
with circular cross-section and continuously varying radius. 

No more than three branches could converge at any point in the 
net_work. This model, however, does not seem to have gained 
the favour of the research workers in the field. 

Fatt replaced each pore with a cylinder in his two 
dimensional network with different grid patterns. He could 
incorporate random assignment of parameters, such as pore-size 
distribution and number of cylinders meeting at a node point, 
in the model. Fluid flow in the network could be followed as 
one pore after another was desaturated. 

All of Fatt's calculations were performed assuming 
a continuous wotting phase in the tubes. As a consequence, 
no wetting phase could be trapped in a desaturation process. 
Dodd and Kiel extended this work by applying desaturation 
steps such that wetting fluid could be trapped. Cylindrical 
pores were assumed to contain only one fluid at a time; thus 
allowing the displaced fluid to be trapped if no continuous 
path to the effluent end was available. Thus at any time 
the pore could be either full of wetting phase or non-wetting 
phase. 



22 

Sin glial modified Patt's model by providing for the 

possibilities of different flow regimes in different channels 
during desaturation process. The four basic flow regimes viz, 
single phase flow, displacement, annular and slug flow were 
allowed to occur any where in the model whenever suitable 
conditions existed. In addition, he substituted triangular 
capillaries for circular shaped tubes tc account for funicular 
flow regime over extended range of saturation. A trapping 
factor, evaluated from the comparison of the model with the 
experimental results, was used to effect trapping of fluid. 
This model has duplicated the relative permeability behaviour 
of porous media reasonably well and we propose to use the same 
model for approximating its capillary pressure vs, saturation 
behaviour. 

Other Model : 

23 

Haring and Greenkom have proposed a parametric 
statistical model, from the theoretical considerations, to 
calculate the macroscopic properties such as saturation, 
permeability and dispersion coefficient in porous media. 

Porous medium is approximated by a large number of randomly 
oriented straight, cylindrical pores with randomly varying 
radius and length. The relationship between saturation and 
P contains two constants which are the parameters of the 
pore-size distribution function. They maintain that the two 
parameter incomplete Beta function can adequately represent 
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most of the pore-size distributions encountered in porous 
media* Non-uniformity of a medium can be incorporated by 
adjusting the magnitude of the two parameters in the distri- 
bution function. Of all the models proposed to-date this 
appears to be the most promising for analytical description 
of flow phenomena as well as mass transfer through porous 
media. 

24 

Pay takes et al, while developing a model for the packed 
bed, point out some in-consistency in Haring and Greenkorn’s model. 
Whereas in Haring’s model the minimum tube size is implicitly 
assumed as zero, the smallest pore-size in a packed bed is 
definitely not zero. As a'consequence, the average pore diaiiwVe-/ 
calculated from this model are smaller than the smallest pore 
calculated from sphere pack model. The calculated pore length 
is about equal to the half the average grain diameter of the 
pack. However, this minor inconsistency in the model can be 
removed by substituting for the finite value of the smallest 
pore-size for packed bed and thereby slightly modifying the 
model. 

In their own model, Paytakes et al assume convergent- 
divergent circular tubes as the basic units of which a packed 
bed is assumed to comprise of. These basic units do not have 
any lateral connection thus allowing fluid flow only in the 
vertical direction. By introducing convergent-divergent flow 
channels, they seem to have accounted for inertia effects in 
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flow path. However, this could not he a realistic approach 
for porous media where interconnections are quite important 
in considering the flow behaviour. 

2.2 Pore-Size Distribution ; 

It is obvious that pore-size distribution is very 

important in describing fluid-flow through porous media. 

25 

Ritter and Drake proposed a method of determining the pore- 

size distribution from mercury porosimetry data. The basis 

of this method is the concept that mercury is forced into 

gradually smaller pores against capillary forces as the pressure 

is increased. The fact, that such information does not directly 

give the pore-size distribution is obvious. However, pore-size 

distribution obtained from Nitrogen adsorption agrees well with 

that calculated from mercury porosimeter data on the assumption 

2 . & 

of parallel tube model. Mercury porosimetry will not detect 
the dimensions of those pores that are accessible only through 
necks narrower than the pore itself. The volume, of these pores 
of limited accessibility will be erroneously assigned the 
dimensions of the largest neck leading to these pores. They 
analyzed a number of natural substances and concluded that the 
pore size distribution can be represented by three parameter 
modified Maxwellian distribution function. Haring’s two para- 
meter distribution function is simpler and should be preferred. 
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Meyer attempted to correct this distortion in the 
pore-size distribution due to the presence of the pores of 
limited accessibility. On the assumption that the pore sizes 
follow Poisson distribution, he calculated the probability of 
a pore of radius r Q not being connected to the pores of radius 
r p,,, r Q . This probability was used to correct the mercury 
porosimeter values for the 'Ink-bottle' effect. The details 
of the derivation and the computational method is given in 
Appendix B. 

27 

DulTi.en proposed to determine ’two dimensional’ pore- 
size distribution by injecting Wood's metal into the sample 
and using micrographic analysis. Apart from the experimental 
scheme, he has not published his results yet. 


2 . 3 Characterization of Porous Media on the Basis of Capillary 
Pressure Curve ; 

The earliest attempt in this direction was made by 

pQ pQ , 

Leverett^ * y who put forward a semi-empirical relation 


J (cl 


(I)* 


( 2 . 6 ) 


~w' ~ 'y'cos 9 f 

where J(S ) = capillary pressure function and is dimension- 

less function of wetting phase saturation 
on the basis of dimensional analysis. Leverett's data showed 
that a plot of J(S W ) vs S w (wetting phase saturation) yielded 
a unique curve which described the capillary retention of wetting 
liquid in clean, unconsolidated sand. 


15 


30 

Rose and Bruce extended Kozeny -Carman treatment to 
relate capillary pressure with porosity, permeability and 
textural constant known as Kozeny constant. The expression 
suggested iss 


( K yh = 

yt 03 o K f J 




( 2 . 7 ) 


where P, 


T 


threshold pressure i.e. lim P 

Sw* 1 


•T 


But these expressions, apart from classifying the porous 
media and in some cases predicting flow behaviour, did not 
give any insight in the distribution of the pore-sizes. 

? "Z 

Recently Haring and Greenhorn have derived an expression for 
wetting phase saturation from P c curve as 


1 

T? * 

rC / 1 \ oC+2 

\ v p* > 


S. 


o^ 


■P* J 

c 


(1 - p 1 ,-)' 0 a ( j ») 
, c c 


w 


f (k)** 2 (1- k )pa( f *> 

c 


(2.8) 


y X P* 

o c 


where P* 


fb 


P 

r T 

threshold capillary pressure 
parameters of pore size distribution. 


Thus in case we are able to estimate and jb from P Q vs S w 
curve, we really have an expression for pore-size distribution 


P(x), 
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x* 

j x* ( 1 —x dx 

F(x*) = ^ — ■ — (2.9) 

J x**' (l-x)^ dx 

o 

where x* = r/r taax , the ratio of pore-size to the largest 
pore size. 

'This. F( x) is then used to calculate the wetting fluid 
saturation, S w , for any P c by the Sq . (2.8). 

2.4 Plan of Attack s 

Prom the analysis of the past work, it is obvious that 

the relation between the po re -structure, its distribution and 

the capillary characteristics of porous media is not well under- 
2 

stood. Meyer assumed Poisson distribution and suggested a 

method for applying correction for ’ink-bottle’ effect. Haring 
23 

and Greenkorn suggested a statistical model for capillary 

pressure which could be useful in mass transfer calculations. 

But they approximated the pores by straight cylindrical tubes 

beside assuming unimodal distribution of the pores. In real 

porous media these are not strictly valid assumptions. In the 

interpretation of P_ curve for pore-size distribution, Patt 

used an arbitrary factor ( {&) which could take the value -1. 

There seems to be considerable disagreement between different 

2 

workers regarding the value of this factor, j' . Meyer had 
implicitly assumed the value 1 for ^ in his calculation of 
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of probability function for the pores of limited accessibility. 

To find a method which can interpret P curves satis- 
factorily, we will do the fo llowing ; 

a) networks will be constructed for different pore-size 
distributions and by allowing them to imbibe nonwetting fluid 
under successively higher pressure, P c curves will be generated 
for different distributions. 

b) An attempt will be made to workout a procedure for 
interpreting P in terms of the pore-size distribution. To 
arrive at the number distribution of the ' po re-sizes from P. 

O 

i 

curve, Patt 1 s model with different values of q/. will be tried 
and the result willbe compared with the actual input to the 
network. 

c) Existence of 'ink-bottle* effect in a real porous 

2 

system will be studied by the method proposed by Meyer . The 
corrected pore-size distribution (volumetric) will be matched 
with the one estimated from the P curve to determine whether 
the effect is significant in natural system. Corrected pore- 
size distribution from the network P c curve will be compared with 
the actual distribution to test the reliability of the approach. 
To this end, different expressions for characteristic volume, 
in Meyer's method, in terms of power to pore radius, will be 
tried and compared with the actual input to the network. 

d) To obtain the pore-size distribution from an actual 

2 3 

P curve, Haring and G-reenkorn model will be used to estimate 
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the parameters of distribution by non-linear least square 
method. Network model will be used to determine how close 
the estimated distribution is to the actual distribution. 


***** 



CHAPTER 3 


MORELS USED IN THE PRESENT WORK 
3.1 Generation of P„ Curve- from Known Pore-Size Distribution-; 

"j. m-mmmmmm. ** mmmmmmmmrnmmmmmmm Q, r ■ > 'ur i i »wii»rri i - ■ r ■ n 11 i;-- i r-mrr r"i: . r u r "'~ ,T '" L r '" Tr 1 ^ 

Use of network model for generation of P curve consists 
of two steps; 

a) Construction of a -network simulating the porous 

material with a known pore-size distribution. This idea was 

originally suggested by Patt and the methodology followed here 

22 

was developed by Singhal . In the network f pores could merge 
at and branchout from a node point. A network of the size 
6x8 was used in our simulation in which the nodes were arranged 

in square grid pattern. Thus at each node point a maximum-- of 

SI 

eight tubes cou-ld meet. Patt suggested, from experimental 

observations that each tube was connected to 10 other tubes 

on. Ike. 

implying that ^ 6 tubes could meet at a point. This fact 

was incorporated in the model by modifying the cumulative 
pore-size distribution obtained from capillary pressure curves, 
so that 25$ of the tubes have zero radius. This was achieved 
by compressing the real distribution between 0.25 to 1.0 instead 
of 0.0 to 1.0 by the following expression; 

1’ = 0.751 + 0.25 

1’ = new value of the distribution 
1 = original value 


where 



20 


Whereas the tube radius and location were assigned in 

random manner by Monte-Carlo technique, the tube length could 

assume only two values viz. equal to the side length of the 

square or its diagonal length. The grid length was arbitrarily 

chosen. The list of the input to a typical network simulation 

program is given in the Table 1. The details of simulation 

22 

has been discussed by Singhal and for the sake of completion, 

the relevant points are included in the Appendix A, 

b) Generation of Ca;;ill’ary Pressure Curves Capillary 

pressure characteristics of the network is evaluated by modifying 

22 

and adapting the computer program written by Singhal Some 
additional checks have been introduced in the program to verify 
that material balance for each fluid as well as for the total 
fluid is satisfied at each stage of computational cycle. 

Initially the network is assumed to be full of wetting 
fluid and. is allowed to be displaced by non-wetting fluid. A 
fixed pressure drop ( z&P) is applied across the network and 
it is allowed to desaturste the wetting fluid in a quasi- 
static manner. This process continues till outflow from the 
network in a particular time step is negligibly small; implying 
equilibrium fluid distribution in the model. Theoretically 
in a quasi-static process the time step should be as small as 
possible. But this would mean larger computation time. Hence 
we used the time steps of the order of 200 /'.sec. initially 
and later increased it as the desaturation process proceeded. 
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Instead of waiting indefinitely for zero outflow A in a time 
interval, we cut short the calculation for a /y P whenever 
the outflow in a time step is less than or equal to of the 
outflow in the first time step. The saturation of the wetting 
phase in the model is then calculated thus obtaining one point 
0 n the vs S„, curve. Pressure drop ( AP) across the model 
is then successively increased and in each case the equilibrium 
wetting phase saturation is obtained. In this way the whole 
P vs S curve can be obtained over any range of pressure and 
saturation. 

Singhal^ used a trapping factor ( 0 , 2/pi^ , where 
X-A- and are the viscosities of non-wetting and wetting 
fluid respectively, to trap the fluids in the pores during 
desaturation. This was merely a correction factor to match 
the model behaviour with actual samples. In static process 
like capillary desaturation, the viscosity of the fluid is 
not of any significance as far as equilibrium saturation is 
concerned. Equilibrium saturation of the wetting fluid depends 
solely on the curvature of the remaining fluid interfaces in 
various tubes. In our simulation of the network the tubes have 
assumed the size equivalent to the volumetric average radius 
of the pores. In real porous systems ( granular) the pores 
are convergent-divergent in shape. As a result the neck radius 
in any tube is definitely smaller than the volumetric average 
radius. As the capillary pressure of any tube is determined 
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by the size of the narrowest constriction in the tube, the • 
actual capillary pressure will always be larger than that 
calculated from the average radius. 

In the present work, the trapping factor is the ratio 
of the volumetric average pore radius to the neck radius for 
each tube, thus implicitly simulating the convergent - divergent 
shape of the pore channels. At every stage the pore-entry 
capillary pressure is modified by this factor to calculate Prp 
at the neck. This pressure is then used to determine whether 
the interface will move through the pore or will be blocked. 

3.2 Correction for "Ink-Bo t tie” Effect ; 

In natural porous systems pores are surrounded by 

2 

smaller pores; recognized as "ink-bottle" effect by Meyer . 

This results in assigning larger volume to the smaller pores 
while the volumetric share of larger pores is less than the 
actual volume. Evidently the use of bundle of tubes model to 
calculate pore size distribution is likely to be erroneous. 

p 

Meyer attempted to correct the measured data for this effect 
by calculating the probability of a larger pore not being 
connected through equal or larger pores to the mercury source 
in mercury poro simetry . He assumed that pore-sices are 
randomly distributed and follow Poisson's distribution. In 
deriving the probability < j Q (r Q ,r) dr, he overlooks the fact 
that the pores are mutually exclusive and asserts that this 
will not introduce any serious error in calculation since only 
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the’ zero’ term in Poisson’s distribution is used for calculation. 
Also, assumed randomness of pore sizes exclude, from consideration 
clusters of very large or very small pores such as fractured 
or vuggy system. 

Other simplifying assumptions made by Meyer are; 

(i) that the centre of an r pore lying in volume v(r ,r) will 
ensure its being connected to an r Q pore, (ii) that only those 
pores which arc isolated from all equal or larger pores will 
fail to fill and (iii) volume of a pore may be expressed as 
Kr" where K is a constant. 

The first assumption is valid for spherical pores and 
seems intuitively sound for others. The implication of the 
second assumption is that isolated groups of mutually connected 
pores are not being considered. This, of course, is a serious 
simplification but necessary for handling the difficult subject 
of aggregates in a random distribution. It is hoped that the 
correction applied for individual isolated pores is the major 
correction needed. The third assumptions implies that larger 
pores have larger length, associated with them. This sounds 
reasonable as the larger pores are associated with larger grains 
for uniform packing. This assumption has one snag that K is 
not a constant. However, K as a constant appears only under 
the integral sign and cancels out. The assumption, in fact, 
is not that K is, a constant but that its mean value in one 
integral is equal to that in the other. 
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Hence, lie proposed that measured pore-size distribution 
should be corrected by this method to obtain actual distri- 
bution for granular porous system. Details of the method is 
given in Appendix B. 

Computer program, written for this scheme, is given in 
Appendix D. 

3.3 Haring and Greenhorn Model : 

p'5 

Haring and Grecnkom ' have suggested the use of two 
parameter incomplete Beta function as the universal distribution 
function for the pore-size in -a porous body. It was claimed 
by them, and rightly so, that almost any unimoaal distribution, 
whether symmetric or skew, could be approximated very closely 
by the adjustment of the two parameters. Normal distribution 
of the pore-size was not favoured for the following two reasons; 
i) the distribution is symmetrical, 
ii) the independent variable varies from -qO to + cp ; 
none of which is strictly true for a porous media. 

One more implicit assumption in their approximation 
is that the ratio of minimum to maximum pore size is zero, 
which is not always true. That seems to he the reason, why 
for a granular pack the average pore size computed on the 
basis of this distribution function comes to even less than 
the minimum pore-size computed from geometrical considerations 
of granular packing. Also, the average pore length is calculated 
to be equal to half the average particle diameter; vfnile ideally 
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it should be equal. This fact has been pointed out by Paytakes, 
.24 

Tien and Turian . However, this inconsistency may be over 

come by slightly modifying the distribution function to; 

r* * iQ> 

\ (x-3^7 ' (l-x) ,v dx 

g(r*) = Y — — (5.1) 

!(x-x ( 1-x) ^ dx 

'o 131 

where r* = dimensionless radius of the pore. The ratio 

of pore radius to the largest pore radius, 
x^ - The ratio of smallest to the largest pore 
radii 

<4 , ft = Parameters of distribution. 

Based on the assumption that the actual pore-size 

distribution can be approximated reasonably by incomplete Beta 

25 

function, Haring et al arrived at the following expressions 
for the wetting fluid saturation, 8 W : 


where 


JL 

(<*+£> + 3)1 

^ = f * ( ..y..+2 ) T /a : 

5 


x^ +2 (1-x)' 8 dx 


= B 


(JL 

V p* 

c 


oC+2 B 



(5.2) 


(5.5) 


B( p* I ci+2 , ji> ) is incomplete Bota function with 
. c 

.y^+2 and jt as parameters.. 

P* Dimensionless capillary pressure, P c /P^ 

Prp Threshold capillary pressure i.e. P c when 


% 1 • 
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In this work we shall try to fit the Eq.(3.3) "by 
non-linear least square method. Standard- IBS program SDA- 
3094 will be used for this purpose. 


***--* 



CHAPTER 4 


RESULTS Mp_ DISCUSSION 

As already mentioned in 'earlier chapters network 
seems to be a good model for predicting the flow behaviour 
of porous media, A part from its closer and better physical 
similarity with porous media, capillary pressure curves 
obtained from it look quite realistic. However, there may 
be some weaknesses in this model which will be discussed 
later in this chapter. 

Capillary pressure (P ) vs wetting fluid saturation (S w ) 

from the networks are given in Figs. 2a to 2f . Each network 

was generated from a specified pore-sizo distribution, which 

was assumed to represent a hypothetical porous medium. Pore- 

size distribution used to generate networks were; arbitrary 

distribution, Beta distribution and Poisson distribution. 

Besides, a couple of bimodal arbitrary distributions were also 

used. The networks with bimodal distribution were used to find 

out if the P curve could give some indication of the multi- 

modality in the distribution. In each case, except the first 

network, the largest pore-size used was 20 JX . The input 

distribution was compressed in each case between 0.25 to 1.0 

22 

to accommodate Eatt's Beta factor as suggested by Singhal . 

Inspite of the very distinct and extreme forms of the 
distributions used to generate the networks, th« calculated 
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P c curves were not significantly different from each other. 

Thus giving the impression that capillary pressure is not a 
very strong function of the pore-size distribution in the net- 
work. Also threshold pressures (P^) obtained from the network 
P c curves (corresponding to S w = 1.0) do not always agree with 
that predicted from the bundle of tubes model. Threshold 
pressure, from the definition, depends on the size of the 
largest tube in the input face of the network. But, in view 
of the fact that tubes are assigned randomly at different node 
points, it is possible that the largest size tube may not 
occur at the input face., In that case I'm predicted from the 
bundle of tubes model will be different from that evaluated 
from the network. Prom the inspection of Pig. I.e, it is 
obvious that the largest tube size (20 /J- ) does not appear in 
the input face and hence the calculated P^ from the network 
is higher than that predicted for 20 Li tube radius ( 0.259 psi). 
Also, during the computation of P curve from networks, more 

V 

than one fluid/fluid interfaces wore observed in a continuous 

flow channel. This factor has not been accounted for in the 

bundle of tubes model and can be one of the reasons resulting 

in the difference in the two values of P^. 

Other features of the network likely to have influenced 

the P curves are? (i) number of tubes (ii) ratio of the tube 
c 

length to the radius (iii) the extent of inter-connection 
between the tubes and (iv) existence of different flow regimes 
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in flow channels. 

31 

(i) Eatt estimated number of pores in a l/2 n cube sample 
of sandstone to bo of the order of 10,000. In our networks, 
the number of tubes varied between 89 to 101. Whether such 
a small number of tubes (pores) could reasonably approximate 
the flow behaviour of real porous media with huge number of 
pores, is open to doubt. Y/ith such a small number of tubes, 
it is the relative position of the tubes which matters 

more thorn the size distribution. If the Investigation could 
bo repeated with larger sizes of networks (500 to 1000 tubes), 
the result may throw some light over thi,f question. Our 
limitation in computer memory and time prevented us from 

A 

trying bigger networks. Eatt asserts that tubes more than 
200 in number do not substantially effect the nature of ? c 
curves from networks. However, he aid not compare any 
experimental P_ curve with that of a network to prove his 
point that tube number is not a critical factor. 

(ii) The second factor, which seems significant, is the 
ratio of the pore length to the pore radius. In our networks, 
this ratio varies from shout 6 to 60 depending on the size 
and the location of the tubes relative to the grid. Earlier 
investigators suggest that the length of a pore varies accord- 
ing to some power of r.. This factor varies between 0 and 1 , 
according to Dallavalle^, Eatt”* puts this value at -1. In 
one of our calculations grid size was reduced to about l/5th 
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of the original value -thus reducing the maximum value' of 
length to the radius ratio from 60 to 12. The two P curves 

v 

(Fig.lh) are not significantly different. Further investi- 
gation with different values of the grid-size is desirable 
before* drawing any conclusion. 

(iii) The third factor to be considered is the Beta factor 

1 ?o 

suggested by Fatt and used by Singhal in his simulation 

of networks. This factor essentially represents the extent 

of interconnection between the pores in porous media. Applied 

to networks, this implies average number of tubes meeting at 

a node point, Fatt , from his experimental observations, points 

out that in a sandstone rock each pore is connected to 10 other 

pores thereby implying tjaat on the average 6 tubes should meet 

at a node point. In Fig* 2h the three P curves have been 

c 

obtained with different values of the Beta factor. It is 
obvious that the two P curves, for Beta factors 6 4 8, are 
not very different from each other, the P curve for Beta factor 

O 

equal to 4 is significantly different from the other two. Hence 
it may be inferred that the extent of interconnection in a 
network, represented by Beta factor, is quite important and 
is likely to be different for different media. In our simulation 
this factor has been taken equal to 6. 

( iv) There is experimental evidence to support the existence 
of different flow regimes in flow paths of physical models of 
porous media where one fluid displaces the other. In such 
systems, the interplay of viscous and capillary forces comes 
into existence and the flow regimes depend on A P and flow 
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rate. This phenomenon is accounted for in the network hut 
not in the bundle of tubes model. This could be another weak- 
ness of the bundle of tubes model-. However, for mercury- 
vacuum system this situation, of different flow regimes in 
flow paths, is not likely to occur. 

Simulated networks used for generation of P„ curves were 

p 

tested for ^oodviess of fUrby "V test.. The pore-size distributions 

of the simulated networks were not found to be similar to 

the input distributions within a reasonable confidence limit. 

We suspect that the scheme of pseudo random numbers used in the 

simulation might not have been truly random. Consequently it 

was decided to use the actual distributions of the networks 

rather than the input distributions which they were supposed 

to represent. Any further work in this direction should ensure 

2 

that the simulated network does pass 'X test to represent 
the input distribution. However, the weaknesses mentioned above 
may not be very serious in view of the results obtained from 
the network model. 

Since our networks have discrete pore-sizes, ideally 
some step type P 0 curves would have been expected. But from 
the limited number of the points on the P c vs S y/ plot, a smooth 
curve could be drawn. It appears that the pore-sizes may not 
be the exclusive factor determining the value of capillary 
pressure in a porous medium. Other factors like connectivity 
etc might also be important. 
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The idea of relating P Q to the size of the pores has 
its root in the bundle of tubes model of porous media. Even 
the model proposed by Haring and Greenkorn^ are the modifi- 
cation of the same model, providing random orientation of the tu 
tubes and consequent inter-connections. Therefore to correct 
for inter-connection between the pores, some modification is 
called for in the interpretation of P_ curves based on this 

• V 

model. 

The effect of "ink-bottle" phenomena seems to be signi- 
.ficant in the networks, resulting in the cumulative pore-size 
distribution (volumetric) obtained from P curves, being 
quite different form the actual distributions. Actual distri- 
butions were computed by counting the number of tubes of * 
different size and calculating their volumes. The difference 
in the two distributions can be clearly seen in the Pig. 3a to 
3- f» Attempts were- made to correct the distributions (obtained 

p 

from P c curves) by Meyer's method for "ink-bottle" effect. 

In using this method, slight modification was made in the scheme 
given in the Appendix B. As is obvious, in a network, tube 
lengths are independent of radius. As a result the pore volume 
as well as the characteristic volume defined in Meyer's method 
becomes proportional to the radius squared rather than radius 
cubed. The corrected distributions (volumetric) is compared 
with the actual and that obtained from P c curves in Pigs. 3a 
to 3f. 
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Effect of the pores of limited accessibility ("ink- 
bottle" effect) on the volumetric distribution appears quite- 
pronounced in some of the natural substances investigated by 
us. The substances included porous plate, fritted glass, pelleted 
gel etc,. Meyer’s method as given in the Appendix B was used 
for these calculations and the result is given in Pig. 4a to 4di ; 

25 

P c data for these substances were obtained from Ritter and Drake 
and are reproduced in the Table 3. The same method was applied to 
P c data of some sandstone samples obtained from literature but 
no appreciable difference in the experimental and corrected 
distributions was observed ; implying that this phenomenon is not 
very significant in those samples. In an attempt to closely 
approximate the actual pore-size distribution, the exponent of 
r in Meyer's method was given a number of values (3 to 1.5). The 
curves with-the exponent equal to 1.5 were observed to be better 

than others in approximating the actual curves (Pigs. 5a to 5f). 

0 ^ 

Physically this means that the pore length is<£r • . This 
distortion from the expected relation of loCr 0 may be due to a 
number of factors, like some bias in our networks, and to a 
greater degree on the extent of inter connection between the pores. 
Attempts were also made to estimate number pore-size 
distribution from T n vs S„ T curves obtained from the networks. 

In the absence of an analytical expression incorporating net- 
work concept of porous media, investigations were carried out 
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to find out if Faffs model, of bundle of tubes, could be 
useful. In that model (Eq. 1.3) Eatt" 1 used an exponent which 

related the pore length to the radius. The same equation, with 
the different values of oC > was used to calculate pore-size 
distributions (number) from P curves of the networks. It is 

* V 

apparent from Pigs. 5a to 5f that for the values of lying 

between -2.5 to -3.5 (mostly -3.0) the calculated distributions 

are reasonably close to the actual distribution. The physical 

— 3 0 

implication of this value- is , that the pore length (l) o<. r 
In our networks the tube lengths were independent of the radius. 
But because of the fact that we have distinctly two groups of 
lengths, side lengths and diagonals of the grid, in our net- 
works, pore lengths could not be said entirely independent of 
the radius. Also, the small sample size, extent of inter- 
connections and the non-randomness of the pseudo random numbers 
used in simulation, might have introduced some fictitious 
relation between 1 and r. 

23 

Haring and G-reenkorn model relating pore-size distri- 
bution (number) to capillary pressure and saturation, was 
tested by estimating the parameters of the distribution from 

P„ vs S , curves by non-linear least square technique. These 
c w 

parameters were put back in incomplete Beta distribution 
function to calculate the number distribution. These calculated 
pore-size distributions did not match the actual distribution 
obtained by counting the number of the tubes from the networks. 
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One of the curves fitted to the network P curve is given 
in Fig. 6 # 


irom the limited number of analysis performed by us 9 
it is doubtful whether this model could at all relate the 


pore-size distributions to P c curves of real porous 

irorn our study it can be said that networks 

used as a valid model of porous media and pore— size 

could be obtained from P^ vs S curves. 

c w 


media, 
could be 
distributions 


*** 



CHAPTER 5 


CONCLUSIONS AND RECOMMENDATION 

The- classical model (bundle of tubes) of porous 
media have a number of weaknesses. As a result the inter- 
pretation of P curve based on this model is likely to be 
erroneous. Network is a better physical model of the porous 
media and P curves obtained from it is quite realistic. 

v 

True pore-size distributions (volumetric and number) can be be‘B«v 
estimated from P curves by modified Meyer's and modified 
Patt’s method ivvcRoodecL dHvs wotk. 

Future work in this direction should ensure that the 
actual pore size distribution (number) of the simulated networks 
reasonably match the input distribution. More work is needed 
to find out the effect of tube length-radius ratio, the extent 
of interconnection, and number of tubes in the network on 
the P curve generated. Experimental work with physical networks 
arc necessary to verify these theoretical findings. 






Fig. la - A typical network used for generation of capillary pressure curves. 
Input data table 1 . Numbers indicate tube rad. in microns. 
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Fig. 1c -Network from arbitrary bimodal distribution. Table 2. 











Fig. 1 e - Network used to obtain P c vs. 5 W curve ( Fig.2e) indicating a b sense of the 
largest size . tubes in the input end. input incomplete beta distribution 
with parameters t-0 J-0 : Numbers indicate tube radius in micron. 





F ig. 1 f - Network for poisson distri 





$w (Wetting phase saturation) 
a~P c vs. $w curve obtained from network I Fig. la 
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D c vs. 5 W curve obtained from network II Fig. lb 
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0 
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6 
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0:27 


9.5 

0.362 


10 
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17 
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*** 
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TABLE 2 

RBI TRA P. Y DISTRIBUTIO N USED IK SIMULATION OF 

NETWORKS 


ore-Size 

Cum. Number 


Distribution 

0.0 

0.0 . 

2.0 

0.0434 

4.0 

0.1304 

6.0 

0.2608 

8.0 

0.3478 

10.0 

0.3912 

12.0 

0.4346 

14.0 

0.4780 

1 6. 0 

0.5650 

18.0 

0.9130 

20.0 

1.0 

0.0 ' 

0.0 

' 3.0 

0.075 

,5.0 

0.325 

8.0 

0.400 

10.0 

0.450 

12.0 

0.500 

14.0 

0.550 

1 6.0 

0.650 

18.0 

0.950 

20.0 

1.0 


#*-**■ 
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TABLE 4 


COT 1 '. ALISON OF EXPERIMENT a 

1 Ah PORE -SIZE DlSTEIBTTTTnTVT 

0 B Y AIT ED FROM MERCURY PC 

>R0 SIMETRY TO TEAT 

CORRECTED 

BY MEYER' 

S METHOD 


Poore !> o_rou_p_ ..Plate 



i "o re - ! :> 1 £5 o ( M 1 c ro n ) 

Mr outer than equal 
to. 

Cum. Fractional. 
Pore Volume 
Exptl. 

Cum. Fractional 
Pore Volume 
Corrected 

1 . 066 

0.00578 

0.09797 

0 , 533 

0.01734 

0.03578 

0.213 

0.53757 

0.71379 

0. 1 066 

0.87861 

0.92197 

0.053 

0.93642 

0,93954 

s — 

CM 

O 

> 

0.99422 

0.99422 

0.01060 

1.00000 

1.00000 

P.yrex UP Fritted Glass 



I’o re-Siae (Mi e ro n ) 

0 reader than equal 
to 

Cura. Fractional 
Pore Volume 

Exptl. 

Cum. Fractional 
Pore Volume 
Corrected 

1 .0 66 

0.00592 

0.10029 

0.533 

0.07101 

0.24686 

0.213 

0.94083 

0.94083 

0.1066 

0.99408 

0.99408 

0.05331 

1 .00000 

1.00000 

0.02132 

1.00000 

1.00000 . 

0.01060 

1.00000 

1.00000 
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Table /] (contd), 
c ) Activa te d Clay : 


Po re - St ze ( Micron ) 
Greeter than equal 
to 

Cum. Fractional 
Pore Volume 
Erptl. 

Cum. Fractional 
Pore Volume 
Corrected 

1.066 

0.00513 

0.08693 

0.333 

0.01026 

0.01604 

0.213 

0.4615 

0.14971 

0.1066 

0.36410 

©.44945 

0.05331 

0.63077 

0.63897 

0.02132 

0.85128 

0.85128 

0.01 060 

Pelleted £ 0 1 

1 . 00000 

1.00000 

Po r e ~ S i z 0 ( M i c ron) 
Greater than Equal 
to 

Cum. Fractional 
Pore Volume 
Exptl. 

Cum. Fractional 
Pore Volume 

Co rrected 

1.066 

0.1685 

0. 16701 

0. 533 

0.22472 

0.51417 

0.21:52 

0.60674 

0.75752 

0,1066 

0.74719 

0.79563 

0.05231 

0.85393 

0.85803 

0,02132 

0.94382 

0.94382 

0.01060 

1.00000 

1 . 00000 
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TABLE 6 


CALCULATION OF NUMBER CUM. DISTRIBUTION FROM P LATA 
OE NETWQB.ES BY LATTVS MODEL WITH MODIFIED EXPONENT 


Network 

P 

V 

s 

Computed Number Distribution . 

c 


G w 

P-.= -2.5 

. = -3.0 o'. 

.=-3.5 

“Til 

L2l 

—I 5-1 

(4) 

(5) 

(TS 

(7) 

3 

0.37 

14.00 

1.0 

1.0 

1.0 

1.0 


0.50 _ 

10.5 

0.6920 

0.5247 

0.3854 

0.33 


0.75 

6.9 

0.5948 

0.3962 

0.26236 

0.201 


1.00 

• 5.16 

0.5155 

0.3086 

0.1 883 

0.116 


1.50 

3.46 

0.3685 

0.1716 

0.0904 

0.047 


2.0 

2.60 

0.2603 

0.0310 

0.011 

0.0044 

II 

0.324 

16.0 

1.0 

1 .0 

1.0 

1.0 


0.5 

10.5 

0.91 

0.831 

0.752 

0.667 


0.75 

6.9 

0.729 

0.558 

0.424 

0.313 


1.0 

5.16 

0.596 

0.390 

0.258 

0.161 


1.5 

3.46 

0.526 ' 

0.315 

°* 195 

0.112 


2.0 

2.6 

0.318 

0,130 

0.064 

0.029 


3.0 

1.73 

0.235 

0.067 . 

0.027 

0.010 

III 

0.287 

18.00 

1.0 


1.0 

1 .0 


0.323 

16.00 

0.70 


0.485 

0.560 


0.370 

14.00 

0, 666 


. 0.43 

0.510 


0.430 

12.00 

0.666 


0.43 

0.510 


0.516 

10.00 

0.538 


0.287 

0.358 
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liable 6(contd) 


(2) 

(3) 

(4) 

(5) V 

(6) 

0. 646 

8.00 

0.354 


0.1184 

0 . 862 

6.00 

0.283 


0.0678 

1 .290 

4.00 

0.221 


0.0365 

2.580 

2.00 

0.1015 


. 0.0 

0.259, 

20.00 

1 .0 

, 1.0 

1.0 

0.30 

IT. 4 

0.8657 

. 0.76 

. 0.72 

0.37 

14.0 

0.7195 

. 0.54 

. 0.44 

0.5 

10.5 

0.6296 

0.418 

0.305 

0.75 

6.9 

0.5143 

0.286 

0.180 

1.0 

5.16 

, 0.431 

0.204 

.0.118 

1.5 

3.46 

0.265 

0.07 

0.030 

2.0 

2.6 

0.260 

0.067 

0.029 

3.0 

1.73 

0.194 

0.027 

0.010 

0.37 

14.0 

1.0 

1 .0 

1.0 

0.50 

10.5 

0.5186 

0.215 

0.30 

0.75 

6.9 

0.5081 

0.202 

0.29 

1.00 

5.16 ' 

0.4032 

0.116 

0.182 

1.50 

3.46 

0.2763 

0.035 

0.07 

2.0 

2.6 

0.2206 

0.0128 

0.03 

3.0 

1.73 

0.190 

0.0043 

0.0115 


( 1 ) 


(7) 


IV 


V 


0.16 

0.090 

0.042 

0.0 
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Table 6 (Contd.) 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

0.259 

20.0 

1.0 

1.0 

1.0 

0.30 

16.0 

0.988 

0.98 

0.977 

0.37 

14.0 

0.8718 

0.795 

0.746 

0.50 

10.5 

0.6216 

0.4428 

0.360 

0.75 

6.9 

0.4009 

0.1818 

0.1144 

1 .00 

5.16 

0.3703 

0.1514 

0.091 

1.50 

3.46 

0.3029 

0.095 

0*054 

2.0 

2.6 

0.2077 

0.0275 

0.017 

3.0 

1.73 

0.1653 

0.002 

0*0056 
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APPENDIX A 

NETWORK MODELS PLOW COMPUTATION 


Pressure Computation s 

Assuming that flow in the pore channels of a porous 
system is essentially viscous in nature (Re 4.75), Poiseuille’s 
law of viscous flow could be applied for flow computation in 
each tube of the network. For single phase flow in a cylindrical 
tube of length 1^ end radius r^, the expression for flow rate, 
q.j_ is given ass 


q i 

or 


7T r i 

~SJZ 


AP 

1 i 


(A.1) 


*i = K 1 

where is the pressure differential causing flow in the 

ith tube. For multiphase flow and h,Pj_ may be modified 
depending upon the relative proportion of the different fluids 
present and the nature of the flow regime respectively. Knowing 
total fluid conductivity, K^ of each tube and maintaining total 
volumetric balance a Kirchoff type equation can be written for 


each node point of the networks 

I k b.L ( + P A = ° 

s=1 ,8 
k=i-1,i+1 


(A. 2, 
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Knowing the values of K g , ? c and the pressure at the 
boundaries, the set of simultaneous algebraic equations can 
be solved by Gauss-elimination technique. 

Can ill ary Pressure : 

Capillary pressure, P Q , is applicable only to those pores 
which undergo displacement flow. For. single phase, slug or 
annular flow regimes, capillary pressure is essentially zero. 

In this computation P c is taken as equal to the pressure required 
to squeeze an fluid/fluid interface from a mode point opening 
(radius r nQ ) to a tube of radius r t ). For a tube of 

triangular cross-section 

B 0 = 1.78 Y ooa e (U P-) (A.3) 

r no r i 

where 'y 7 = interfacial tension and © is the contact angle which 
in our case is the receding contact angle. 

Plow Behaviour in a Single Tube s 
i) Single Phase Plow: 

The total volumetric flow rate through a tube of 
triangular (equilateral) cross section, with the radius of the 
inscribed circle, r, viscosity of the fluid /U , and pressure 
gradient AP/ is 

q _ 1 ( A. 4) 

jr ^ * r 

ii) Annular Plow: 

We used modified Yuster relationship as suggested 
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22 

by Singhal for flow computation. 

A 


q. 


•w 


% = 


1.35 

• sP- 

1.35 


r 

-Aw 

si' 

An 


AP c2 

“ w 


ZM 




A 


n 


2S t~ + 
n 


(A. 5) 


f£(1-2 


A 


£))!■ 

*- 

W / 


(A.6) 


wher© S denotes saturation and subscripts w and n denote wetting and 
and. non-wetting phase respectively. 


(iii) Slug Plow; 

Total flow rate is given by 


q 


1.35 


tot v r 3 K A + + 6a/ V 


ZiP 


v > V* vv xx XI w / ^ 

where a is the side of the triangular cross section. 

( iv) Displacement Plow; 


^ = 1 t (A ' 8) 


Selection of Tube for Plow Change ; 

Depending on the availability of fluids at a node point, 
transition between flow regimes will occur. In many instances, 
channels leaving a node point may be capable of carrying more 
of a particular phase (fluid) than is available at that time. 

In such cases, saturation and flow regime changes in the tubes 
will have to occur to unload all available fluids. A table 
of possible transition between flow regimes is given in Table A-1 . 
Some times there may exist more than one pore which might be 
undergoing flow in the regime to 4 be altered. Selection of tubes 
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for flow changes in such cases is done by Monte-Carlo method 
with a bias for flow capacities of the competing tubes. 

General Scheme for flow Computation ; 

Words in block capital letters are subroutine name 
referred to in the program listing. 

The network is generated in the program by SHU with 
specified pore-size distribution and grid size. This subroutine 
also computes the volume of the pores connected to each node 
point and flow conductivities to the two fluids. • Other data; 
fluid viscosities, interfacial tension, pressure differential 

m across the model and contact angle, are read in the main 

** ' . ■■■ . : - ; : ' ; , ■ ' , ; ' . ' 

program CALC. , The radius of the opening at the node points, r n , 

is chosen to be slightly larger than the largest tube radius 
of the network. 

The next step is to determine the pressure distribution 
in the model. Pressure at each node point is computed by 
solving the system of equation A-2 by PRSOL through SOLVE. 

Apart from updating the conductivities of each tube after every 
time step, SOLVE also immobilizes those tubes which may cause 
outflow from the first row and inflow from the last row. Also, 
those tubes, whose capillary pressure is greater than the* 
pressure differential available and are oppositely directed, 
are blocked by the action of trapping factor. When a tube is 
blocked its conductivity is set equal to zero and the system 
of equations A-2 is solved. This is repeated till no more 
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tubes are blocked. While calculating P for tubes in the 

c _ 

first row, the node point opening is assumed to be infinite. 

For other rows the specified value of r n is used for the 
calculation. 

The total outflow of each phase (fluid) from each node 

point is calculated under the existing pressure distribution 
* 

by FLOW, Outflow from a tube is entered as inflow to the 

connected node point by CQED(l). The imbalance between inflow 

and outflow of a particular phase at a node point is corrected 

by CHANGE which also effects the flow changes necessary for 
* . . ' : : 

this correction. Tho capillary pressure and flow conductivities 

for the corrected flow distribution is computed by C0ED(2) and 

stored. Outflow from the model is obtained by adding together 

the inflow to the last row of node points. The volume of the 

fluid left in the model is computed to evaluate the saturation. 

This process is repeated till the equilibrium wetting phase 

saturation in model is obtained for that pressure drop. 

Steps of Computation ; 

(i) Call SIMU to simulate the network model. Assign 
tube sizes, locations etc. Calculate the volume of each tube, 

. its conductivity to displacing and displaced fluids. 

(ii) Call SOLVE to update the conductivities through 
C0ND(2) and solve the system of equations A. 2 by PRSOL to 
evaluate the absolute pressure at each node point. In case 
some of the tubes are blocked,’ call CQNL(l) to zero the 
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conductivity of those tubes and again call PRSOL to solve 
the systemof equation A-2. Repeat till no more tubes are 
blocked. 

(iii) Choose node points in descending order of absolute 
pressure and call FLOW to calculate outflow of each fluid from' 
each tube and node point. 

(iv) Calculate imbalance at each node point to test 
whether overall material balance of incoming and outgoing fluids 
is satisfied. 

(v) Imbalance in the material balance of each phase is 
* •• removed by assigning the excess incoming phase to various 

tubes by Monte Carlo method through CHANGE without forcing flow 
change in the tubes. 

(vi) Assign* outflow from each tube as the inflow- to the 

connected down —stream node points* 

(vii) In case of back flow, calculate inflow to those 
tubes and node points in FLOW. If there is any imbalance either 
in total fluid at the connected node points or imbalance in any 
phase in any tube, call CHANGE to assign the differential volume 

to various tubes. 

(viii) Calculate outflow from the model by adding the 
inflow to each node points in the last row of the network. 

(ix) Calculate total fluid volume left in the model by 

i -£* J-Vv t y\ p pi cii tnlo© at tli© end of tiie 

adding the volume q£ the fluid in eacn 

time interval. 
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(se) Calculate the saturation for each phase. 

(xi) Repeat steps (i) to (x) till there is practically 
no change in the total outflow from the model in the successive 
time steps (less than Yjo of the outflow in the first time step) 
for the same C> P. 

(xif.) Test whether .the sum of the cumulative outflow and 
the volume of fluid remaining in the network for each fluid 
at the end of every time step is constant. 

(xiii) Repeat the sequence of operations with higher AP, 
to get the value of equilibrium S w (wetting phase saturation) for 
different P c . Thus P c curve is generated. 
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TABLE A-1 

TRANSITION BTOEEN BLOW BEffTTyms 


Plow Regimes 

To / Prom 

Single 

Phase 

Displacement 

Slug 

Annular 

Single Phase 

X 

1 * 

1 * 

1 

Displacement 

2 

X 


- 

Slug 

- 

- 

X 

1 ; 

Annular 

- 

-]** 

1 

1 


x No change in flow regimes. 

1 Easy transition, 

2 Transition influenced by capillary forces. 

Transition not possible without going through 
other flow regimes. 

* Modified single phase flow for non-wetting flow 
because of a film of wetting fluid at the pore 
walls. 

** Possible only' for non-wetting phase displacing the 
wetting phase. 
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APPENDIX B 
MEYER' S METHOD 


The radius r of a pore is defined to be equal to the 
radius of a capillary tube which will just alio w mercury to 
enter it against capillary force at the same pressure as that 
at which the pore will fill, if connected to a mercury source. 

Let }J ( £ ) be the fraction (volumetric) of the total 
void space occupied by r-pores- i.e, pores whose radii lie in 
the range r to r+dr. The. number of such r pores in a unit 
volume will be given by • 




« * v ■. 

a(r) dr = — MXi- 5 * (B. 1) 

Kr p ' 

where f = fractional porosity 

■ac 

Kr = vol, of a single pore of radius r 
It is necessary to designate one representative point 
in the interior of each pore which we shall term "centre of 
the pore". Bor our purpose, any point will do so long as it 
is not too close to the wall of the pore and is consistently 
chosen, e.g. eg of the void space will serve. Also, we define 
a volume dv- so small that the probability of the centres of 


two or more r-pores being located in this volume may be 

neglected. .1 vpl'A X ■■ ■ 

Let the probability of m r-pores being located 
to any volume v be tf> m (v,r)dr, Then for randomly distributed pores 
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( ' v+dv ’ r ) dr = HiJj(v,r)dr ft -(adr) dvl 
U^_ 1 (v,r)dr (adr)dv 


(B.2) 


which on re-arranging gives 


d ( K fm dr ) \S a dr ) T J _ ' , 

<!v = 5! exp[_-(adr) v] (B*3) 

Now we define a volume v(r Q ,r) such that, if an l-pore 
is located in v(r Q ,r), it may be expected to be connected to 
an r 0 -pore. This can be done if wo assume that two pores are 
connected if the centre of the smaller lies in a volume 
adjacent to the larger pore, which is the difference between 
the volume, of a pore whose radium •'equals the sum of the radii 
of the two pores minus the volume of the larger pore i.e, for 

r o> r 

v(r Q ,r) = K(r 0 +r) 5 - Kr Q 3 = K(3r 2 r Q + 5r r 2 + r 3 ) (B.4) 

Now if in Eq.(B~3) v is replaced by v(r Q ,r), we get the_ 
probability of m r-pores all being joined to the same r Q -pore, 
The probable number of such groups of pores which would exist 
in a unit volume would be given by the product of. this pro- 
bability times the number of v (r Q ,r) volumes in a unit cube 
(which is equal to the number of r Q pores in a unit cube). 
This number is equal to 


a(r )dr (v(r 


, xv , r y aj 


a(r Q )dr (^ 0 >^) 
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We are here interested in ^(r Q ,r)dr which is- the 
probability of an r^ pore not being connected to an r-pore. 

Let a porous rock be subjected to mercury pressure such 
that all pores of radius ^ r would fill if they had access to 
mercury, and let f (r) be the measured fraction of the pore 
volume that is filled up with 'injected mercury at that pressure. 
Let ©(r) bo the actual fraction of the pore-volume devoted to 
pores of radius ^>,r. 

If the simplifying assumption is made that only pores 


of radius ^ r which will fail to fill are those that are not 


connected to any other pores of equal or larger radius. That 
is the probability of an r pore failing to fill is given by 

O * - or 

the multiplicative law of probability ^ (r ,r)dr. 

; r , ' 0 ° ... 

The total number of pores of radius^ r Q which will fail 
to fill is given hy the expression 



Hence \ (r) 


= f a(r)dr 

sj' 

CO 

[j kVo (r o- r) ar ) 

(B.5) 

r o 

r 


r 

■' r (r. ,r) -< 


= 0(r) Li - 

* 0 ! 

ZJ& ■ J 

(B.6) 

J a(r)dr 



This is the expression from which we shall obtain true pore- 
size distribution ©(r) from the actual distribution f (r). 


Numerical Approximation ; 

Since the above expression is not convenient for numerical 
calculation the following approximation are used. 'The radius r 



85 


is discritized so that it takes only a finite number of 

values r n / r n _ 1 /r., and 

Z.\r i 

a i = 

r. 


Ji-1 " r i 

j~' 1 a(r)dr 


1,2 


~ 1 ^ r i^ “ . \ ^ r i^i ^ 
v i = e(r i ) - 0(r i-i ) 

'tlj = ex p l- a i v ( r i’ r j)] 

.t( r i» r j ) = °k T\ fkl 

k =1 1=1 


*§■ 


Substituting these expressions inEq.(B. 6) gives the approxi- 


, . £ TT v Ha 

— - . JL _. r - k=i 1=1 i {v 

i' ~ ^Yk = 2L.^k x L 1 T~ ] (B ' 6a} 

X ®k 

k=1 


mate equation 

n (X- *1“ : ■***” V *} 

r n - ) = "T Vi_ = ) h. x i 1 - K ~ 1 

k=1 k=1 


Algorithm : 

(1) Assume a value of and calculate a i and _Yij* 

(2) Using these values solve Eq. (B.6a) for )j^. 

(3) If computed value is not close enough to the 
as sum bdy value, mkke another assumption and repeat 

d ■ ■ ■ : d : v ' ' ■ V ■ ■ 

this proofs. 
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API EMU IX C 

GLOSSARY OP TEMS TTSPP 

Porosity (f) is the ratio between the void space and bulk 
volume in a porous medium. 

Capillary Pressure (P Q ) .is the difference in pressures across 
a fluid/fluid interface. Capillary pressure curve describes 
the changes in capillary (threshold) pressure of a parous 
as a function of its fluid saturation. 

Thro sho Id Pro s sure (P^) is the minimum excess pressure must 
be maintained in the non-wetting phase to enable it enter 
a given porous medium. 

Saturation ( S ) of any fluid in a porous medium is the ratio 
between the volume of the fluid in the medium to its pore 
volume. 

Contact Angle (Q) is the angle (through the displacing phase) 
that the fluid/fluid interface makes with the surface of the 
porous medium. 

Wetting Phase is the fluid in a porous medium whose contact 
angle with the solid surface is less than 90°. If the contact 
angle is more than 90° it is called non-wetting phase. 


**** 
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U 

SIGMA 

THETA 

POR 

R(I) 

PC 

SGAMA 

ESTHU 

HU 

DIOTU 
SV(I, J) 
PSI(I, J) 
PSIR(I,J » ■ 


APPENDIX D 

PROGRAM LISTING OP MEYERS METHOD 

Humber of data in the set 

Interfacial tension 

Contact angle of non-wetting phase 

Porosity 

Pore-size corresponding to P(I) 

Capillary pressure 

Cumulative Volumeof Mercury Injected. 

Assumed value of nu 
Computed nu 

Difference between estimated and computed nu 
v(r- .,£•)■» characteristic volume 



ISN 


SOURCE STATEMENT 


FORTRAN SOURCE LIST 


ss 


0 SIBFTC MAIN 

C MEYERS METHOD 

C PC IN PSI 

C CAT*'. MEANS THE USUAL VALUES OF CHANGE ARE NOT TO BE USED 

1 CONNONMi/PSI (20,20 3, TEMP (20,201, SUMAt 20 ), SGAMAt 20 

2 COMMON/ A2 / RS ( 20 ) , RC (2 3, SIGMA 

3 COM M0N/A3/NU (50*3 ) , FOR 

4 DIMENSION P ( 20 ) ,SV ( 20 ,20 ) , R( 20 ) , A ( 20 ) ,0 IFNU ( 500 ) , E STNU ( SCO ) 

5 REAL NU 

6 READ 6 y , N 

lb 60 FORMAT (14) 

11 PRINTX20 

12 120 FORMAT ( 1H1 ) 

13 REAC62, SIGMA, THETA 

14 62 FORMAT (5F10. 3) 

15 65 FORMAT (5F1Q.5) 

16 THETA* THE TA*3. 1416/130. 

17 RE-AD62 , POR 

20 PRIM50, THETA 

21 50 FORMAT (5X»*THETA=*, F10.5) 

22 N=N +1 

23 0012 1 = 2 , N 

24- R( I 3*0.0 - 

25 RS C 1 1-0*0 

26 RC ( I ) =0.0 

27 A( 13=0.0 

30 S U M A ( I ) =0 . 0 

31 00 12 J=1,N v 

32 SV ( I , J ) =0.0 

33 PSI(I,J)=C.O 

34 P S I R =0 . 0 

35 12 CONTINUE 

40 00 101 I K = 2 , N 

41 101 REAC65 ,P( IK) ,SGAMA{ IK) 

43 PRINT20 , (P ( IK ), IK=2 ,N ) , (SGAMA(IK), IK=2,N) 

54 DG 40 I M=2 , N 

55 40 SGAMAt IM)=SGAMA{ IMSJ/PCR 


57 

60 


P R I IS T 5 

5 FORMAT ( 5X» X2HACTUAL VALUE ,10X,2HIK 


, 5X, 10HTRUE VALUE 




0Q1C J=2 , N 

.CAT*OiO..: . ' 

CIFMJ( !)»;=. O . .T ; ; 

I K * 2 « . : . 

ESTNU( IK)=SGAMA< J) 

2 QALL SI(IK,PSIR,A»ESTNU, THETA, N,P,R,SV,J,J) 
IF( IK. EQ.2.AND. J.EO .2)PRINT2v ,(R(I ),I-i,N) 

IF ( IK. EQ. 2. AND. J.EQ.2 (PRINT 20, ( (SV(I, IJ ) ,1 J 

20 FORMAT ( 5X , 32 3.3) 





CGG145 


ISN 


133 


136 


141 


144 


147 


152 


155 


156 


161 


162 

25 

163 


166 


171 


172 


173 


174 

7 

175 

6 

176 

1* 

200 


205 

8 

206 

22 

207 



CGG14 5 


set PCS STATEMENT 


FORTRAN SOURCE LIST MAIN 


I F ( ESTNU ( IK) .GT .O.OC1 ) CMANGE = C.OOI 
IF(ESTiMU(IK) .GT. .0 1) CHANGE=O.Ol 
IF ( ESTNUI IK) .GT.o.l >ChANGE«0.1 
I F ( E STNU ( IK) .GT.0.9 JChANGE=0.C05 

1 F (C I FNIJ ( IK ) .GT .0.0 ) ESTNUI IK + I) =£STNU( IK J + CHANGE 
I F ( C IF NU { I K ) . LT.t ■.!) ) LSTNUI IK+1 )=E$TNU( IK ) -CHANGE 
I K = I K+ 1 

IF ( IK.GE.5C0) GO TO 7 

GO TO 2 

CHANGE=CHANG6/2 . 

1 F ( Cl FlilJ (IK) .GT .0.0 ) ESTNU ( IK + I )=ESTNU( IK J+CHANGE 
I FtCIFNLM IK) .LT .0.0 ) ESTNU ( IK + 1 ) = E STNU { IK J-CHANGE 
I K * I K+ I 
CAT*:." 

GC TO 2 

PRINT 6,SGAMA(J)»IK,eSTNU(IK) 

F0RNAT<5X, FiC.5,UX* I5,5X,FIC.5) 

CONTINUE 

PRINTS* U(J) ,SUMA(J),J = 2,N) 

FCRMAT(I0X*E15.7,10X, £ 15.7) 

STOP 

ENC 

IBMAP ASSEMBLY MAIN 


NO MESSAGES FOR ABOVE ASSEMBLY 


'■* 


CGG145 

ISN 


SGLPCt STATEMENT 


FORTRAN SOURCE LIST 


0 

X 

2 

3 

4 

5 

6 
7 

10 

13 

15 

16 

17 

20 

21 

22 

25 

26 
27 
30 

33 

34 

35 

40 

41 

CGG145 


$1 13F TC SUE 1 

SUBROUTINE VQL ( P,R , THETA, IK, N , SV»M ,IM I 

■ CONKWAl/PSI (20,20 ), TEMP (20,20 ) »SUMA( 20) »SGAMA (20 
CQMM0N/A2/RS (20 ) ,RC (2 ), SIGMA 
CCMN ON / A 3 / N U ( 500 ) , FOR 

DIMENSION P( 20) ,SV(20,2G) ,R(20) , A (20 ) , D IFNU ( 500 J 
F ACTOR =l./fab947.3*lC.**4 
DO 2 1=2, K 

R( I >=FACTCR* 2 .*SIGMA*CGS(THETA)/P{ I ) 

R( I)=ABS(R (I ) ) 

RS( I )=R ( I )**2 
I RC { I ) { I ) 

R ( 1 ) = 2 *i • C 

RS ( 1 ) *R ( 1 ) **2 
N1*N+1 
DC 6 1=2, N 
CO 6 J * ?. t N 

I F(R ( J) .GT.RC I) )G0 TO 7 

5 SV( I,J )*2.*R( I) *R( J)+RS( Jl 

GO TO fa 

7 S V ( 1,J)=2.*RU)*R( JJ + RSU) 

6 CONTINUE 
D01C 1*2, N 
DC1C J *2 , N 

10 SV ( If J)*A8S(SV( I,J) 1 
RETURN 
END 

IBMAP ASSEMBLY SUB1 




NO MESSAGES FOR ABOVE ASSEMBLY 




CGG145 

ISN 


SOURCE STATEMENT 


0 $IBFTC SLE '1 


FORTRAN SOURCE LIST 



IBMAP ASSEMBLY SU82 


CGG145 


NO MESSAGES FOR ABOVE ASSEMBLY 
CC-G145 . ' , ■ 


000000 


W.'k 1 ' 


*** OBJECT 


SUBROUTINE Si ( I K,PS IR ,A,ESTNU,THETA,N,P,R,SV» IL, JK ) 

COM NUN /Ai / PS I <2 , 20 ) * TEMP ( 20 , 20 ) , SUM A ( 20 ) , SGAMA { 20 ) 
COMMON/ A2/RS { 20 ) * RC (20), SIGMA 

CCMMON/AB/NUCSCK 1 ) , FCR 

Cl MENS ION P(2>n ,SV( 20,20) ,R(2C),A( 20 ) , DI FNU ( 500 ) ,ESTMU(5< 

REAL NU 

CALL FL UN (99999) 

N1=N 

MU ( I K ) =ESTMU ( IK } 

SUMA (1 ) = 0 . U 

IF ( JK.EQ.2 )CALLVOL( F,R, THETA, IK,N, SV,M,IM) 

I = I L 

A( I )=NU( IK )*<1./R(I )-:./R(I-l))*PQR 

CO 12 J=2,N1 

B = - A ( I )*SV(I , J) 

GO TO 1 
GO TO 2 


IF (B. GE.88.0) 

I F ( P . L r . (-88.0,) ) 

PS1 ( I , J ) *£>P ( B) 

GO 10 12 

1 CONTINUE 

3 PSI (I, J)-1.0E30 
GO TO 12 

2 CONTINUE 

4 PSI (I , J)-0.0 v 
12 CONTINUE 

J = JK 

DO 8 K * 2 , 1 L 
TEMP ( I , J )»U. 

PRCC = 1 . 

DO 10 1=2, JK 
DELTAR=R(L-1 )-R(L) 
PRCC=PROD*PSI(K,L) 
IF(PROQ.LE.0.1E-37 5 FR0D=0.0 
10 CONTINUE 

8 TEMP (I , J)»TEMPU»J)+A(K)*Pfft)D 
PSIR*TEMP ( I , J ) 

9 CONTINUE 

SUMA (I ) “SUMA ( 1-1 ) +A ( I ) 

5 FORMAT ( 5X, 8E15.7 ) 


RETURN 

END 


l'&"HRS. 21 MTS. 
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AF'MBIX E 

PROGRAM LISTING NETWORK MODEL 

The first subscript in each variable in the program 
refers to the node point. 

X(l) to X(8) The radii of the 8 tubes 
X(9) to X(l6) Length., of the tubes 

X(l7)to X(24) .Fluid content and flow regimes in the tubes 

1.0 wetting (one phase flow) 

2.0 non-wetting (one phase flow) 

3.0 Displacement, (non-wetting close to the 
node point) 

3.5 Displacement (wetting closer to the node 
point) 

4.1 Annular flow (displaced phase wetting) 

4.2 Annular flow (displacing phase wetting) 

4.5 Slug flow 

X(25)to X(32) Flow or no flow situation 

1 .0 Flow occurring 

2.0 Flow blocked 

1.5 Flow provisionally blocked 

2.5 Flow permanently blocked 

X(33)toX(40) Flow conductivities of the tubes to displaced fluid 

X(4l)toX(48) Flow conductivities of the tubes to displacing 

fluid 

X(49)toX(56) Volume of the tubes 

X( 57) toX( 64) Displaced fluid content of the tube 
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X(65) to X(72) 
X( 73) 


X(74) 

X( 75 ) to X( 82 ) 
X(83) 

X(84) 

X( 85) to X(92) 

X( 93) 

X(94) 

X(95) to X( 1 02 ) 

X( 103) to X(110) 

X( 1 1 1 ) to X( 1 18) 

X( 1 19 ) to X( 126) 

IE 

JR 

GAMA 
THETA 
RAH 




Effective pressure drops across the tubes 
Incoming fluids 

1.0 displaced 

2.0 displacing 

3.0 Both 

4.0 none (node point isolated) 

Absolute pressure at node points 

Appropriate conductivities of the tubes 

Volume of incoming displaced fluid 

Volume of incoming displacing fluid 

Effective capillary pressures 

Provisional total outflow of displaced 
fluid 

Provisional total outflow of displacing 
fluid 

Provisional inflow of displaced fluid in 
tubes 

Provisional inflow of displacing fluid 
in tubes 

Provisional outflow of displaced fluid 
from tubes 

Provisional outflow of displacing fluid 
from tubes 

Ho. of rows in network 
Ho. of column in network 
Interfacial tension 
Contact angle 
Tube radius 
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COM 

Cumulative distribution 


GRID 

Grid size 


TRAP 

Value of trapping factor 


HYST 

Hysteresis in contact angle 

(=0) 

TIME 

Time step 


I STAR-1 PIN 

No, of time steps for each 

pressure 

AY 

Displaced fluid content of 

the network 

ICP 

No, of pressure steps 




ST AT LML'JT 


SUURCl 


fqrtran SOURCE LIST 


2 / 


jIBFTC CALC 3 HP - A 

'COHMOM/SYI'lnKG/y ( 5. , 16 ) , AP { 5!) ,8 ) , GC, VO, VW» TI ME , VOF < 5 , 3 ) , I J , I K, JR » JPP , . .v 
1XA CP! I/.' ) jHYSTrCAT, GAMMA, J , FLO f 5 *: » 8 ) , GRID 3P°00o*0 

COMMON / PRE SS / NN , A { 5 G» i>0 ) * AB ( 50) 3PP0-3v3 0 

COMMON / SO L/ A A ! 5 », 5 ’ ) » BB ( 5 J ) , A V ( 40.) , WR ( 4 0 ) , WOT 4 ) ,WI (4 ) , D£-N( A ) 3PPG 5 

COMMON/ SUL B/ CAP { 1 s , U ) , I , I CHANG » AR AT , TRAP I3PP0 7 

C Li M M : ) N / c. i I A N / 0 i F J j W 0 IF 8 PPG 0 0 8 -j 

CUHMUM/CHI N/ ICAn! 50 ,8 ) , AWD I F ( 50, 8 ) , AD I FO ( 50, 6 ) , BASE 3 p^OO j90 

COMMON / S3LC/A03»CP, JIM, I CP BPPv.i * 

COM MOM / S I M/C JM( 15 ) , RAD ( 15 ) *KL ,LK BHPO .. » 

C READ AND PRINT ALL THE RELEVANT INPUT DATA 3PP0Ui20 

CALL FL UN ( 99 9 99) 3 HP (. • 4 . 2 .. 


DC 9 7 IKJ=j.,2 

READ ! , I R, JR, GAMMA, THETA, CP, VO, VW 
PRINTOUT 
93) FORMAT (HU ) 

PRINT 4 , IK, JR, GAMMA, THETA, CP, VO, VW 

READ I U , K. L , L K 


3 P P ■ ' 4.1 
BP PC 
B PPOOi 3 
3 PP<* 13. 
B PPOOi 4 j 
3 PPOOI bv 


PRINT iJ, IP, JR, KL , L K 
READ A- , ( CUM ( I ) , 1 = 1 ,LK) 

PRINT 4'’, ( CUM (I ) , 1 = 1,LK) 

READ 40, ( RAD ( J ) ,J = ). ,LK) 

PRINT 4/ , ( K A D ( J ) , J=i,LK) 

PRINT 4 ,VU,VW,CP 

READ V.i, GRID 

PRINT 40, GRID 

READ i 2 , XA ,HYST , TRA P 

PRINT 12 » XA, HYST , TR AP 

READ 7 , T I Mfc 

PRINT 7, TIME 

READ i. , ISTAR,I FIN . 

PRINT iOtlSTAR, IFIN 
I =1 STA R-l 

READ 3 7 , WR ( 1 + i ) , DEN ( I +1 ) , STW, STO 
I F ( 1ST AR. SO. I ) 1=1 
READ 37 , A V ( 1 ),WO(i) , WO (I), Will) 
A00= WO ( I ) 


7 FORMAT ( IX, F15. i J ) * 

10 FORMAT! 213, 6F10. 5) 

11 FORMAT ( 20X, 4l 4) 

12 FORMAT (8F1<». 5) 

13 FORMAT ( 6F1J. 5 ) 

H FORMAT U 2 F; .0.5) 

Vj FORMAT ( IX, 8F10* 5) 

K FORMAT ( IX, 21 8) 

50 FORMAT { 20X, 15H0I L SATURATION*, F19.4,17HWATER SATURATI 0N=, F2 ). 

51 FORMAT ( 20X»7HTHE PT .» 14,25HIS BETWEEN iAND 4 IN C73 > 

52 FORMAT! 20X,7HTHE PT ., I4»12HHA$ 4 IN C73 ) 

100 FORMAT!//// 8(9X,IH )) 

HI FORMAT ( 3 OX, 21 HP RES S UR E DISTRIBUTION ) 

U7 FORMAT ! //IX, 3U< , 18 H FLUID DISTRIBUTION ) 

31 0 FORMAT ! //20X , I 4,7F13»S// > , 

T59, ,;F OR fi A T IS 0 a » * C U N T'A G'T AN<?UfJ Fo * l , » ill ,, 

?S0 FORMAT ( T. i X , * V I S C 0 S I TY RATIO! OIL/WATcR) **, F£».2) 


B P P 0 0 1 6 
BPP'j .7 
H P o 00 ). o j 
B PPOOI 9C 
3 PRO" 1 7 j . 
BP PC) 2 j. 

B P P 00 2 2 0 
3 PROP 23' 

3 P P 0 J 2 4 
BPPO025 
3 PP00260 
BPPGu27 J 
BPP00280 
BPP00290 
3PPOC30. 
BPP0G31 G 
3 PP 0032 0 
3PP0C-3 3J 
BPP0G34S 
BPP0335? 
3 PP 003 60 
3 PROP 3 7 1 
B PP00380 
B PP00390 
3 PP00400 
3PP094i.> 
BPP00420 
3PP0D43 j 
BP p 0 : 144 ,; 
3PP0045D 
BPP0046U 
3 PPG- ‘47 ■ 
bpp s v: 4b 
BPP0 r '49 . 
3PP005LC 
3 PPL , 6 
3 P P C ' 5 V 


SCJ'JRGh ST A T - M.Nl 


FORTRAN SOURCE LIST CALC 



■l/o'/ \i\ 


7f> 1 F QRM A T ( ’ H .1 ) 

,j=IR*.J R 
JIM=N- JR 

THETA is CONTACT ANGLc AND GAMMA IS INTERFACIAL TENSION 
A RA T = V J / V N 

ARG=TH E T A* . .. t.l o/ 1 6C. 

GC = I • 7 o * G A M 1 1A* C OS ( ARG ) / 68 947 . 3 1 
CAT = 0 . 0 

PRINT 759, TH TA 
READ OAT A FROM tape- 
call SI MU 

PRINT? 7 , ( ( -'( J,< ) ,K=i, d) , J=I,N) 

PRINT3 7, ( (X(J,<) ,K= 57,64 ),J = 1,N) 

DO 2 ' lJ = i,N 
1 K = I J / J K + 

JK=I J-JR*( IK-* ) 

1F( JK.LT.- ))K = IK-1 
I F ( JK.LF. «)JK = JR 
EP( IK, UK) *>:{ IJ, 74) 

CAP ( I K , JK ) = . » 

DO 39 I A* , 0 
X(IJ»IA+16)-,.0 ' 

X{ IJ,I A + 24) = * .0 
X C I J,IA + 136)*HLTA 
1 F ( X ( I J » I A ) * LQ. .)OJX( IJ, IA + 24) = 2. C 
X( IJ, IA+40 )*X(I J, IA+32)*V0/VW 
X ( I J,I A+56)=X( I J , I A +4 3 ) 

X(IJ»IA+74)=X(IJ,IA+32) 

X( IJ, I A + 84 J*.*) 

X( IJ,IA+94)=0.D 
XU J,I A + 102) =0.0 
X ( I J , I A + 1 J, )=.),,) 

x(u,iA+lla)=j. J 

39 CONTINUE 

X(I J, 3 3) =0.0 
XU J,34 ) =0. 

XU J, 93 ) =0.0 
XII J, 73) =1.0 
X(I J, 94) “0.0 
29 CONTINUE 

PRINT I?, GAMMA 
S T W = 0 • 0 
ST03. "J 
.o- A V ( 1 ) = t - . 

""DU 900 I CP=I » 12 
■■ WKlhii.O 
DO 9 TO I CP«1 » 13 
IF( ICP.Gt;.3)IIME = , .001 ' 

I F ( I CP • G il • 5 ) T I M b = 0 . CO 5 0 A 

I F U CP . G L . o ) T I M c = 0 • Cl . 

R AD .’, CP 

PRINT7 , CP , 

DO 700 ! = 1 STAR, IFIN 
IFUCP.CQ.U WOU) =0.0 
IF-{ICP.cQ.UWm) = 5.C 


A P ? T 5 3 
3 PP ’’ 5-i-D 
3 PPO PCiU 
3 P P i't 
8 P P 0 5 7 
8 PP' ; 5dJ 
i PP0J590 
3PP-; 6 . 

HP PC) U 
3 PPO )bcV 
3 P p 0 D 6 s 0 
DPP (. . r>_ j 
t? P D 0 ' ;j 3 2 
3 PP 006* J 
i PP 00 d 5 0 
3 P P 0 o 6 : .• 
BPP', D? I 
o P P 00 6 0 0 
3 P P 006 90 
3 P P n , 7 .« 
8PP‘* '7 s. . 
BP POO 720 
BPP J 1 7 u ■ 
3PP0.74.- 
8 P P OiJ 7 5 J 
B PP00760 
BP PO 77 
8 PPOj79.« 
BP POO 7 90 
B PP00 800 
3PP00810 
BPP 03 8 Z i) 
BPDQ0330 
BPP0J&40 
3 PP 00850 
BP POTS 60 
BPP 00870 
3PP00880 
3PPO’>893 
BPP 0 191 
BPP 00 920 
BPP 00930 
3PPCJ094J 
8PP03950 
BPP 00 960 
BPP 00 97-., 
3PP00980 
BPP0398I 
B PPG0982 
BPP 00 96 3 
BPP0A99, 
B PP Oi 00 J 
3 HP 01 Oi u 
B PP 2 
B P P C .1 1, ■ 3- 




FORTRAN SOURCE LIST CALC 


SOURCE. STATEMENT 
*8 CONTI MUD 

; FRA Mr EXAMINE AND SOL 7c THE SYSTEM OF 
CALL SOLVE 

1 F ( ICP.LT. C) 30 TO 1\j7 
i,u print 

: PRINT PRESS 01 STRIBLT10N 

DO 1 07 IP = , i R 
PRINT * *u 

PRINT i. A- » ( . P ( 13 , JK ) ,JK=*,JR) 

177 CONTINUE 
PRINT x,..7 
A00 = .'» J • 

A W U “ J • . 

SWL= . . i 

I F ( I.Nc.i ) SWL = 1 .0-AVm/AVd ) 

AV ( I ) = v . C 
ILK = j 
S AT L=0 . 0 
SATO=0 . 0 
DO 3.33 KR S = » » N 

3 CALL POINTS IN PROPER SEQUENCE 

IF(KRS.LE.lLK) GO TO . 00 
DO 305 K=KRS,N 
IK=K/JR+I 
JK=K~JR*( IK-1 ) 

IF( JK.LL.0)IK=IK-1 
I F ( JK.LE.O)JK=JR 
IF( JK.EU.JRJGU TO 320 

DPE = EP ( IK » JK !-EP ( IK ,JK+i)+CAP { IK, JK) 
1F(DPE.G£. 0.0)30 TO 320 
305 CONTINUE 
320 I LK=K 

MK= ILK-KRS+ *. 

DO 350 KR = l » MK 
I J = I LK-KR + 1 
J = I J 

IF { ICP . LT. 15 ) GO TO 2 
PRINT 10, I J 
2 CONTINUE 
1 K = I J / J R + „ 

JK«I J- JR* (IK-1) 

IF( JK.LE.O) IK=IK-1 
IF! JK.LE.O) JK=JR 
l:jA=IJ 

• V'T ¥ V' m> 

' ; VWI'-0*O 

VOIDS'*-®.' 

MMK = T , 

BAS E=0 . 0 

IF(X(J,?3) .r'Q.A.O) GO TO 29i 
C COMPUTE FLOW RATES FOR OIL AND WATER 

CALL FLOW . . | t Wl 

VOI =X( J ,83 ) 

VWI=X( J,34) J 

I F ( I K . E Q . i ) VOI = \ . 2 " 


H* P P f f . * 

EQUATIONS i3PP0:u5J 

5 PP 0. j6u 
3PP' 1 

BPPQj - i 
BPPO 0 r f. 

3 PP 01 0 90 
BPPU 1 

BPPu. 1 

6 PPO,,., i! a 
3 PPC, 3 

3 PP 9 1 j ■'+ 
BPP 5. 

3 PP O. 1 & 6u 
3 PPM.-. . I 
3 P P 0 1 i 8 
tiPPOl, 4.9 . 
BPP01200 
3 PP 0.121 
3 P P C jl 2 c 
&PP<U 2 In- 
B PP01 2 4u 
3PP0 25 
B P P C 1 2 6 
BPPO 1270 
BPP01280 
BPPO: 29. 
BPP013.J 
DPP01310 
3 PP 01 32 0 
3PP01330 
BPP0134T 
BPP01350 
3PP01360 
3 PP0137 .< 
(3PP0138) 
B PP01381 
3 P P 0 i 3 9 .»>. 
BPP01391 
BPPO 1 4)0 
BPP01410 
BPPU4L 
BPPO1430 
BPPQ144 i 
BPPO. 45 .. 
BPP0146U 
BP P 0147 j 
8PP01430 
SPP01490 
BPPO" 5 ; . 

INCOMING OUTGOING BPP0151 

3 P P C '■ 520 
3 PP 01 b:> ij 

6 PP.,:54 
BP POT 51 



FORTRAN SOURCE LIST CALC 

SOURCE ST AT;;; MEN t 

1 F ( I K • L U * j. 3 V W 1 = 

1 F ( I K . L Q . 1 3 GO TO 2i4 

UPDATE X (73) =LI QUIDS ..NTERING THE JUNCTION POINT 
I F 1 V 0 1 . c Q » * . AND. V WI . EQ. 0 . ! ) ) X 1 J , 7 3) = 4.3 
IFlVDI.EQ. . . A NO . V W1 . HE . 0 . j ) X ( J , 73 ) =«7 . Q 
IF1V0I.NL. . . AND. VWI .NE.y.olXl J, 733=3.0 
I FI VOl .NE.O.O. AND.VWI . EQ . Q . 0 3 X { J » 7 3 ) = 1 . D 
VGI =A3 $ ( VO I ) 

VWI=ABS( VWJ 3 
IF1IJ.3T.J1 M 3 GO TO 29 
COMPUT . OUT FLOWING VOLUMES 
214 VQQ=X1J,93) 

V W U = X 1 J, 9*- } 

VOO=ABS ( VOO) 

VWO = AB S ( VWO) 

DfcLT4=V0I+VWl-V00-VW0 

1 F ( VO 0 . L r. . jE-i ■.* 3 VO j=D . 3 

IF(VWO.LE.O.: 0E-10) vwo=o.o 

IFIVWl -LE.0.10E-10) VW 1=0.0 

IF1VQI.LE . . IvF-Ij) VO 1 = 0. u 

IF(ABS(DLLTA).LE. C.i — lO) DELTA = .* .3 

DE LQ=V UI -VDO 

DLLGW= VWI-VWO 

DELO=ABS 1 DELO 3 

DELOW* ABS 1 DLL GW 3 

1F( DLLO.LT .Q.1E-1Q) VUG- VOl 

ifidelow.lt. ,ie-i: ) vwo=vwi 

IF1ICP.LT. I 5) GO TO I 
8S8 PRINT37,V0I,VU3,VWI , VWO, DELTA 

I FI DELTA. NE.O.O 3PRINT37, ( X ( J, KA 3 , K A=75 , 82 3 
IF1DELTA.NE.0.0 3PRINT37, (XI J,KA) ,<A=65,72) 

I F{ DELTA. Nr.. . .0 3 PR I NT37 , 1 X ( J , KA + 94 ) , KA= i , 32 3 
IF! DELTA. EQ.C . 3 PRI NT3 7, 1 X 1 J , KA 3 , KA=75 , 32 3 
IF1 DELTA. EQ.O. 3 PRINT37, IX 1 J, KA) , KA=65 , 72 3 
1 CONTI NUL 


2 / 


8PPC, ^ 6 
BPPOi. 570 
3 P P 0 3 PC 
BPP 01 59 
6PP01 5 
d P D C ) 6. J 
3PP03 62 0 
3PPC16-: 1 
BPP 01622 
6 PPQ ! 6',v 
3 P P C i. Or 

BPP -.j: S5 

Q P P 0 .1 6 6 

B PP01 67 0 
3 PPG... 6& . 

BPP C I d 9 
3PP01 7 >0 
BPP 01 710 
8 PPGJ 7 2 
BPP03 73 
BPPOi 740 
BPP01 750 
BPP Ci 7b 
BPP01 77 i 
3PPQI780 
3 PP 01 790 
3PP018.U 
BPP01301 
B PP 01 8 1 0 
3PP01320 
3PP01830 
BPP018V) 
BPP Oi 84i 
3PP01 342 
3 P PCI 65 


UPDATE OIL CONTENT OF TUBES 
DO 21 K=i , 8 
DEN 1 K 3 =X 1 J , K + 55 3 
21 X(J,K+56)=VQF(J*K) 

DO 470 J P = I , 8 
ICAN1 J » JP 3 =0 
AWDI F( J , J P 3 = 0 .0 
ADIFOl J, JP 3=2.3 

IFlXUfJP + 16) .50.0*01 PRINT 11,1 J, JP 

470 CONTINUE 
CASE=0 . 0 

IF ( IK. NE. 1 .AND. VOO. EQ. 0.3. AND. VWO. EQ.Q.03X 1 J, 73) =4.2 
I F ( IK. EQ . I • AND. VOO* EQ . 0* 0* AND. VO I » EQ • Q • I). AND. VWO • E Q • 0 
1 EQ. 0, 03X1 J, 733=4.0 
IF (XI J,73).£Q.4.03 CASE=1.0 
IF(X1J»73J .EQ.4.3) GO TD 270 

COMPENSATE FOR IMBA ANCE IM INFLOW-OUTFLOW AT A POINT 
MONTE CARLO TECHNIQUE 

" IF1IJ.LE.JR) CA$E=! .0 
IF1IJ.LC.JR) GO TO 48 5 


BPP01851 

BPP01852 
3PP02 855 
BPP01357 
BPP0186U 
BPP01870 

Bppoiac 

3PPQ189J 
BP PCI 903 
DPP 05. 910 
3 P POX 92 u 
3PP0193,. 
.0. AND.VWI. BPPC194B 
BPP0I950 
BPP CU 96 
BPP0197 

BYPICKING TU3PP0T 480 

8 PP 01 990 

BPP 07 . ■ 



SOURCE STATEMENT 


FORTRAN SOURCE LIST CALC 



IF (VO I .GS.Vui 
I F ( V j I .CE.Vj: 
1 F ( (V01+VWJ ). 
WDI F = V/ W I -VWJ- 
DI FO = V jl-V.ib- 
1F{ WDIF.Li::., , 
i F ( DIF j.Lt .0. 
I F( CIFj.Lc. . , 
CALL CHANG,. 

1 F ( D I F 0 . IM C . . 


.An u < 

, A N u > 


Vhl .G= .VWOJCAS 
VW1.G2.VW0) GO 
LE.0.0 ) GO TO A BO 5 
■( J.::LT4*VWl)/{ VOI+VWI ) 
(DELTA *VCJI »/ (VOI+VWI ) 
it-. L . ) WDI F = , 0 
) Cl F0=0 * 0 


= 1 . 

TO 


■ u 


-*N 0 * WD I F • L b 


) GO TO 43 


WDI F.NE.Q.O ) PRINTS? » <X{ J,KA+74) ,KA = i,ii) 

TRt AM DATA HAS TO BE CORRECTED DUE TO FLOW 


1005 

l'H>5 


' . 0 R * 

"LANS UPS 

80S 00 1 0 3 K= , » 3 

XCJ,K+56) =r»L.O(<) 

X( J,K+> 6)=FL0( J,K) 

I F { X { J » K + 2 4 ) .GE .1 .5 )GU TO 
IF (X( J,K+6M .LF.C. r )G0 TO * 

X ( J,K+ tiO )=X( J, K + 94 )+X{ J,K + 56)-V0F (J,K) 

X( J,K+i.tC )-)(( J, K + i02)+V0F(JfKJ“X(Jf< + 56) 
X(J»K+56}=V0F(J»K) 

X5 CONTINUE 
27 i CAT =;.* 

CAL L COND(l) ' 

CAT=u. 0 

IF(MMK.FQ.„ . AND. BASE. EQ.C.D) GO TO 291 
I F ( BAS £ • £Q • 1 • u ) GO TO 479 
I C A N ( J * M ) = 0 

GO TO 4B3 

479 DO 481 I C A=6 > b 

IF(BASC.EQ.i.O) PRINT 11 , I J , I C AN ( J , 1 ) , I CAN { J , 7 ) , I C ANf J , 8 ) 

I F ( B4SE.EQ.1 .0) 8ASE=2.0 

M = I CA 

IF( ICA.EQ.6)M=I, 

MMK=I C AN ( J * M ) 

ICAN( J » M ) =U 

IF(HMK.EQ.( ) GO TO 482 
WD I F = A WD I F (JfM) 

DIFO=ADIFO( J,M) 

IF(X( J,M + 64) .GT.0.0 )CASE=2.0 


BPP01. . L 
BPPOi 3 . :■ 
6 PPOI i.‘4 ( 
3 PPG? 

B p p r. i . 3 
3 PP Ot 0 3b 
3 P P 0 i U Im 
3 P P 0 s 
b P P G * . . . 

& P p 0,- ... 0 i 
Cl 13 PPG/ j.; 
b P P G? ~ '< 

B P P C 2 i 2 i 
BPP02.3C 
3 PPG?? 4 
3 PP 02 i 5 
BPPC2*60 
BPP0217O 
BPPG2.B. 
BPPC219 
BPP0220u 
3 P P 02 2 1 0 
BPP0222 ■ 
BP P 0223 I 
BPP02240 
b PP 02 2 5 0 
3 PP 02 cb 
BPP0227' 
B PP 0? 28 0 
3 PP 02290 
3PPO230D 
BP PC?. 31 ') 
BPP02320 
3 PPQ2 33 0 
BPP0234,' 
BPP0235I 
3PP02360 
3 PP 0236 2. 


IF{ CASE. LT. 2.0) J=MMK 

IFCWDIF.LE.O.O. AND. DI FO.LE. 0 . 0 ) GO TO 482 

DO 478 K = „, j E 

FLO( J» K)=X< 4*K+16) 

V G F ( J } K) = X(J»K+56) 

AP(J f K . 

IFCX( J.K+64) . GT.0.0 )AP(J,K)=X{J,K + 64) 

473 CONTINUE • 

430 I F { WDI F. Lt . J . IE -l ) WD IF= 0 . n 
I F ( D I FO • LE .0 . iE -1 3 3 Cl F0 = C . 

I F { D I F 0 • L E • 0 • 0 . AND . WD 1 F» L E • 0 • Q ) GO TO 485 
CALL CHANGE 

ENTER CHANGES IN VOL OUTF LOW/ I NF13 W, AND FLUID CONTENT 
43 5 CONTINUE 
43 6 DO IDS K = 3 > 8 

>U J,K + 56)=DEN(K ) 

X ( 0,K+ l ) = F L 0 ( J » K ) 


3PPC237 ■ 

BPPG233 0 

BPP02.390 

3PP024 

BPP024C 

BPP02420 

3 PPG2430 

3PP0244, 

BPPO2450 

BP P 02460 

BPP02470 

3 PPG? 4o 

BPPG2+9 

B PP 01 50 0 

3 PPO, 5 .. h 

BPP-r’ D, . 

b h p «.• r . s : 



FORTRAN SOURCE LIST CALC 



IF(X(J»K + 24) .GE.1.5 JGO TO 3. '5 
I F ( X ( JiK + 6 * ) . LE . i.O )G( J TO .05 
X( J»K+i i 0 ) =X ( J, K + 94 .) + >' ( J, K+56 )-VOF { J, K) 
X( U,K+..' 8 ) - X ( J,K+„ , 2) +VOF { J , K ) -X ( J , K+56 ) 
X( J,K + 56)=VlJF( J,K) 

CONTINUE 


oPPOl 5+j 
3 PPG ;■ 

3 p p :• -• 

B P P 0 ; 5 7 

a PPG/. tPU 
3PP0: 3 9 

8 P P C 2 6 ■ 
jpp«;'.r-.. 

3 P P 0 s. c L . 

3 P P C > ’ 6 , ' 

3 P PC I 64 
3 P P C 7 5 5 , 

B H° 02 Ot 0 
3 P PC 2 6 7 
BPPC/6R 
B PP 02.69 0 
BPPOZ700 
3 P P (j 2 71 
BPPG27 c .. 
BPP027 yj 
3 PP 02 740 
3 P P 0 2 7 4 j 
1BPP02742 
BPP02743 
3PPG275C 
BPPQ276 
BPPC277 • 
BPP02780 
3PP02790 
3 PP 028 07 
BPP02831 
BPPQ2810 
3PP0281 i 
3PP02820 
BPP02B3G 
B PRO? 8+0 
8PP0285 . 
BP P 02 6 6 
3PP02870 
3 PP 02 86 0 
3PPCV 89- 
BPP029JJ * 
BPP029I 0 
8 PP02 9i i 
3 PP 02 9 13 
BPP0292.," 
6 P P 02 9 0 
3PP02 93 j 
B P P 1 it 94 
BPP0798 
3 P P 0 2 9 o 0 


I F ( MINK . FQ. 0* AND .BASE. i:Q. 0.0 ) GO TO 291 
I F { B AS S • L Q . 1 . 0 ) 33 T3 479 
I C AN { U » M ) = 

4B2 CONTINUE 
481 CONTINUE 
483 BASE=D« 0 
MMK = 

J=I J6 

291 IB=IJ+JR 

DO 333 K = 6 » 9 
IA = K 

1FU.EQ.9) I A = 1 

AVI I ) = AV ( 1 )+X(U, IA+56) 

IF(X{J,IA+56).LT.0.t.UR.X{ J,lA+56) .GT.XI J, IA+46) ) PRINT* J, J , I 4 
I F { X ( J, I A+565.LT.0.€.QR.X( J, IA+56) . GT . X ( J , I A+48 ) ) PR1 NT37 , X ( J , I A+4 
1),XCJ, IA+56) 

330 CONTINUE 

IF{ I.EU.l .AND. I CP. EC. .. ) VQLPOR=AVQ ) 

C COMPUTE TOTAL OUTFLOW OF FLUIDS FROM THE SAMPLE 

IF { VOI .LE. . II E-10) VOI=O.Q 
IF(VWl.Lt.O.iOE-lO) VWi=0.0 
IF(IB.LE.N) GO TO 292 
IF( ICP.LT.i5 JGO TO 3 
PRINT 37 » VOI i VW I 
3 CONTINUE 
AOO=AOO+ VO I 
AWO=AWO+VWI 
DO 290 I A = 7 j 9 
K=I A 

IF(K.60.9) K=*X 

$ATL*5ATt+X( J*K+56) 

satc=sato+-:u,k+48) 

290 CONTINUE 

292 CONTINUE ; 

C PRINT FLOW 

IFIICP.LT. I 

PRINT37, (X(J,KA >,<A=57,64> 

777 PRINTI2* ( X< J »KJ ,K = i7»32I 

IF(X(J,73).E0.4.3) PR I NTS 2* I J 
CONTINUE ’ . " ■ V';v 

I F ( I.EU.l. AND. ICP.EQ.UGO TO 2222 

GU TO 79 
2222 CONTINUE 


SSt! 




FORTRAN 


SOUR 


I F ( X ( J » K + j 5 ) .tQ.i. ) X(J,K+16)=3. j 

212 C0UT1MU! 

213 IF C IK. Me. 2) GO TO 79 
DO 2 5 I A = 7,9 

K= I A 

I F ( K • E U . 9 ) K=_ 

I F ( X { J , K + i. 6 ) . l Q . 1 . 0 ) X ( J , K + l 6 ) =3 . 5 
216 CONTINUE 
79 CONTINUE 
;:( j,-o ) =o . o 
X. ( J , 6 A ) = 0 . 0 
X(J,93)=; . 

X { J , 9A 5 = . 

VO I = , . i 
VWI =0. 0 
29b CONTINUE 
350 CONTINUE 
330 CONTINUE 

PRINT37, AVC I ) 

SO = A V ( I ) / VOLPOR 
SW=1 . D-SQ 
PRINT 50 » SO* SW 
690 CONTINUE 
WI { I)=AWO 
WI { 1) = WQ(1)*ARAT 
CUML= ( STO+ST W ) / AV ( l ) 

WO(I ) = AOO 
STW=STW+WI ( 1 5 
STO = ST O+WO { I ) 

PRINT 37, W0(1) ,W0< I) , ST3, STW » WI { I ) 

W=W0(1 ) *1) . V I 

I F { I . S T . 2. • AND. A 00 . L £ . W ) GO TO SCO 
730 CONTINUE 
800 CONTI NUU 
930 CONTINUE 

CALCULATE AND PRINT POROSITY AND PERMEABILITY 
A JR= JR 

AIR=IR 

PQR=AVIl)/{ ( AJR-1. ) *( AIR-1 . ) *GRI D**3 ) 
PERM=WI(U*VW*( AIR-1. )/( TIME* GRID*! AJR-1. ) *CP ) 

PRINT 12*' POR, PERM 


B P p i( 3 - : 

-3PP03"' 

3 P P 0 ; , 

BP PC? 

BPPOr Re : 

BPP 05 05 
3 p P C ^ 6 

B P P 0 3 - 7 
B P P D :> I. ■ 

b PP 03 

3PPC3 * .. • : 

R P P C > i i. 

B P P 0 3 * 2 • 
bppo:-;.. 1/ : 

8 PPOrU b-L 
3PPC- ; c 
BPP03i60 ; 
3PP0317 j 
3PPU316 ■ 1 

BPP0M9. : 
6PP0320Q ; 

3 PP 032 10 
8PP0322 . 

R P P 0 3 2 3 . i 
8 PP 032*0 
3PP0325. 1 

3 PP0326 ,■ ! 

BPP0327? 
BPP032C0 
3PPG329,. ; 

BPP03303 
8PP0331 I 
BPP0332U J 
3PP0333. 
BPP08340 ; 

BPP 03350 
8 P P 03 3 o 0 
BPP0337 . | 

BPPQ3 58 
BPP 033 90 
BPPG3A0U ' 
BPP‘j3R-„ ! 
BPP03420 
1 2/D +i 


IBMAP ASSEMBLY CALC 




ERROR 2 31.. I S i-i 


|Q VARIABLE REDEFINED IN ITS JWN RAN 


/ERITY W 


._/ / 


SOURCE STATEMENT 


FORTRAN SOURCE LIST UB I 


CALL PDIF 

AVOID INFLOW FROM OUTLET END OR FROM INLET END 
UBES 

K=I J 


DO 76 K L = 2 j 5 

I F { K.GT.JIM) GO TO 76 
COL=D. 

I JK=I J 

IF (X(I JK,KL+24) .Gt. 1.5) GO TO 75 
AS=0. j 

I F ( K.GT .INN. AND. K.LE .JIM) AS = 1."> 

IF I AS. £0.1.0. AND. X( lJK,KL+&4) .LT. 0.05 AS = 2. 0 
I F( AS. EQ.2.U.AND.KL .G 7.3. AND . KL . LE . 5 ) C0L=0 . 5 
IF ( K.LE . JR. AND. X { I JK, KL+64 ) . L T • • J ) COL = 0.5 
1 F ( X { I JK» KL + 1. 6) .Lfc.2.0) GO TO 55 
IF{X(IJK,KL+16) . G T .4.0)00 TO 55 
IFIK.Lu.JR) GO TO 55 
D P = 0 . ) 

DP=X(I JK t KL+84) 

DP=DP* ( x.+TRAP) 

A D P = X ( IJK, KL + 6 4 )-X ( I JK, KL+84) 
IFIDP.LE.U.G.AND.ADP.L&.O.O) GO TO 55 
IF( AOP.GE.O.O.AND.DP.Gt.O.O) GO TO 55 
DP = AES ( DP ) 

ADP= AS S ( ADR ) 

I F < ADP .LT.DP ) ClJL = : .5 
55 IFtCOL.EQ.O. 5) X { I J K, KL+24) =1 . 5 
I F ( COL . EQ. 0 . 5 ) C 0 = 0 . 5 
76 CONTINUE 

IF (C3.EQ .U.D ) Gl) TO 77 
CAT =0. 0 


COUfiT = COUNT + 1 

J = I J 

CALL GONDI I ) 
77 CONTINUE 
B 7 CONTINUE 

DO 82 I L = 1 , 8 


I F ( X ( I J ,73) . E 0 • 4 . G ) GO TO 82 
If (XU J,IL+24).GE.l .5) GO TO 82. 
CPA=X( IJ, IL+84) 

CPT+XI U,lL + 74) *CPA 
VS= i+XUJ*IL + 74) 

JK)=CPA 



DO 

XUJ,i. . . . . — 

:•:{ IJ,<+202)=0.0 
X ( I J , < + 1 1 0 ) = 0.0 
XU J,K + 118)=3.| 

I A= IK ' , -? 1 

J a = J < 

IF(K.EQ.l) i A= I K-l 


■i P P C •> , 1 

TRAP APPROPRIAT, T3 PP •' ,.VJ 7 

HpP'J h ■ 

A PPC401 1 
3PP04U/ 

3 P P ;j : + z, , ; 

b ? P \ 1 *• 4 
BPPQU5 • 
3PPC4 6 
3 P P 4 1 
A P P •. ’■* 

5 P P O'f I- Ju 
GPP 

6 P P G + a. i, 

3 P D QU LA 
3 P P04L : 

3 P P \ih 4 
&PP<>4x5 
3PP0*18U 
3 PP 

3PPG417 
BPPG4U 
3 PP 04i 9 ■ 

3 PP04 2 j 
3PP042i 
BPP0422v. ‘ 

B PP 042 3 0 j 
3PPG425. 
BPPQ426 ’ 
BPP04270 
BPP0428Q 
5PP043U 
BPPC431 i 
BPP 0432 D 
3 PP 0433 0 
3PP0434 » 
BPPQ435 
BPP0436U 
3PP0457 
BPPQ438D 
BPPQ439 j 
BPP 04400 
BPP044i, 
BPP0442U 
BPP 0443 . 
BPP 04440 
3 PP 0445 < 
BPP044o : 

B P P 0447 0 
BPP 04 4 SO 
GPP 0449 
DPP. 

3 P P 0 
3 P P ' 



FORTRAN SOURCE 

LIST UB .t, 



SOURCE ST 4 IF Me NT 






1 F ( K . L t . 3 ) J 6= J K + 1 




bpp 


IFIK.G. .5 .AnO.K .LE. 7) JA=JK-l 




cap -tS; 


1 F ( K . 35 • 3 . AN 0. K . Lfc . 5 ) 1 4 = I K + 1 




BP pc ■ 5 6 ' 


IFCK.Sc.7) 1 A=1 K- i 




?■ ? P ■„ 5 ? 


I JC= ( I A-l ) * JR+ J A 




BPPO-,5 •/. 


A A { 1 J j I JC ) =-X { I J,K + 74 ) 




f P P 0 > ' ■> 'J 


IF(IJC.LE.JR) V A = YA +AA ( I J, 1 JC )*CP 



3 P P v- A 6 

7 

CONTI MU : 




BPP A.-o . 


BB (IJ) = - VA-C PI 




3 P P C 46 : 

B 5 

4A ( I J, 1 J ) - VS 




3PP0C63C 

3 6 

CONTINUE 




3 p P 6'. A 6- 4 


I F ( IC.N".). .AND. COUNT. Q.J) GO 

TO 

78 


B P P v 46 j 


DU 66 IN = I,'1N 




B P w 0* ,: 66 . 


INR=1N+ JR 




BPPQV, M* 


AB( IN) =BB( 3 NR ) 




BPP C ~ o -• 


DO 66 JN= 1 » NN 




b P P . 1 6 9 


JNR=JN+JR 




RPPOV/uU 


A ( IN, JN)=AA( INR, JNR ) 




0 PP047 i, o 

6 6 

CONTINUE 




BPP 047 . 


CALL PRSOLt i ) 




BPP047 ) 


CALL PRSOLC;) 




hPP 0- 7.- j 


DO 75 K = 1 , N 




B PP 0475 0 


8 B ( K ) = C P 




3 PP 0476 j 


lF(K.GT.JIM) B3(K)=U. 




GPPC477 


IF(K.LE.JR.OR.K.GT.JIM) GU TO 

73 



BPP 047 SO 


KK=K- J R 




BPPQ479Q 


BB{ K ) =A3 ( KK ) 




BPPC48, .. 

73 

IK=K/JR+l 




BPPC481 ‘ 


JK=K-J R* (IK-1 ) 




B PP048/:0 


IF(JK.cQ.O) I K= I K- 1 




BPP0483.-;. 


IF( JK.EQ.O) JK= JR 




BPPQ484G 


EP(IK, JK)=RB(K) 




3PP0485:* 

75 

CONTINUE 




BPP04860 

78 

CONTINUE 




3PP0437 


PRINT it), COUNT 




BPP0488 J 


I F { COUNT . NE . J ) GO TO 80 




BPP0489J 


IFUC.EQ.i) GO TO 80 




BPP049UC 


RETURN 




BPP049IL 


END , 




BPP 049 L I 


U- \ . IB MAP 

ASSEMBLY 

UB I 

3. 2 / 0 4 

' S'i. 






f *'■ ; 






-ES 

FOfe", ABOVE ASSEMBLY 








SOURCE STATEMENT 


FORTRAN SOURCE LIST 


it'- r/ 


SIRFTCSUB 3 

SUoR iJTIMw. Flow 

C 0 M M 'J N / S Y M A R G / X ( 5 0 , I b j ) , A P ( 5 0 , 8 ) , & C , V 0 , V W , T I M E » V 3 F { 5 0 
». XA *l.PU. , i ” ) t HY ST, CAT, GAMMA, J , FLO { 5 , B ) , GR1 D 

C FLOW CAP I CITY OF TUBE S 

ARAT = VU/VW 


N = I R* J R 
IB=N-2*JR 

X { 0 t 9 b ) = * 

X ( J » 9 4 ) = j • 

UU 200 K=i « 3 

V0F(J»K)=X(J»K+56) 

mP ( J, < )=’/. 

I F { X ( J j K + 64) .GT.0.0 )AP{ J,K)=X( J t K + 64> 

X ( J > K +• 9 4 ) = 0 • 0 

X( J,K+L' 2 )=■'•./ ' 

X(J,KH10)=i.i . 

X(JtK+ fc ]8)*0.0 

IF(X( J t 73 ) .FU.<+. ) GO TO 2. ,0 
C CASE 2. M>. AMS BACK FLOW 

CASE = 3 . ■*' 

F LO { J t K ) =X ( J » K+ j. 6 } 

IF(K.Gt.3.AND.K.LE. 5) CASE =1.0 

IF (CASE.EQ. I.A. AND. AP ( J t K 3 . EQ . O.D > CAS£ = 2.D 

IF(CASE.EQ.Z.u) AP( J,K)=X( J,K + 64) 

IF( X{ J ,K + 24) .GE .2.0 )GU TO 200 
VA = AP( J»K) 4X(Jj K + 7 4) 

VA=VA*TIME 
VA = ABS ( VA ) 

IF ( VA. EQ.O.Q )G0 TO 200 
I F (CASE . EU .2 .* ) GO TO 5 

GO TO 7 

5 I F ( X ( J, K + 16> .EQ.4.2 ) FL0<J,K)=4.2 
I F { X { J » K+l 6 ) .EQ.4.1 } FLO ( J » K ) =4. I 


mm 


IF(X( J.K + 16) .EQ.3.3 )FLO( J»K)=3. 5 
IF(X( J,K+16) .EQ.3.5 )FLO( J,K)=3.0 
7 IF(FL3{ J,K).LE.2.0) GO TO 10 

I F( FL3( J,K ).LT . 4.0. AND. FLO ( J, K) . GE .3. 0) GO TO 20 

IFCFL3(J»K).GE.4.3) GO TO 30 
10 IFIFL3 ( J»K } . EQ. i . 0 } XI J,K + 94)=VA 
IF| FC3( J,K) .EQ.2.0) X ( J , K+I02 )=VA 

GO TO 9 ; 

2 IF(FL3( J»K).,1Q.3.5) X< J»K+94)=VA 
I F ( FLO { J » K ) . EQ. 3 . 3 ) X ( J » K + 1 0 2 5 = V A 

30,' GO T3 40 

I F CJRL3 ( 3# K ) . EQi 4-i ) X ( J * K + 94 ) =VA 
IF(FL3(J*K).dQ.if.2) X ( J * K + 1,02 ) = VA 
S0= X( J j K + 56 ) / X { ‘J » K + 43 ) 

n>W""X. *” o lJ 

AKRW=X(J» K+4S SW/ARAT+S W#SW*( 1.-2.Q/ARAT ) ) 

4KRN=X( J,K+3Z) *S3*SQ 
AAK = AKRW/VW+AKRN/Vtihi' . 

AKRW=4KRW*X{ J,K+74> /AAK/VW 
AKRN=AKRN*X( J»<+74l /AAK/V0 


bp p < : :.? 

d P P 0 5 : - > 

} , I J, I R » J K i 3 P P G : . 

3 F P 05 ->• . 

B P P i ■ 5 : 

3 P P G > ■+ . . 

3 pp r 3 t*. , 

8 P F C5 4 3 
3 P ° -Li ; 4 r, 
l"i P P 0 - +5 i. 

3 P 0 v 2 4 o . 
LPPOb-T 
Bpo 05-''s, j 
3 P P 05 4 9 0 
3PP033 
UP P C5 .3 . 
BPP0552G 
3 P P 05 5 5 0 
3 P P C 5 5 -V 
BP PC 5 5:: 

3 P D 0 5 5 6 -j 

3 PP 055 7 0 

3PPC5 5 3 

BPP0359 

3 PP 05600 

BPP056I- 

6PPQ562 

BPP0563 

B PP05640 

BPPG565. 

BPP05660 

BPPC567 > 

3PP05680 

3PPG569*. 

BPPC57 o >.* 

BPP05710 

3PP05720 

BPP05730 

BP P 05 7 40 

BPP05750 

3 PP 05760 

BPP05773 

BPPQ57L 

8PP 05790 

8 PP 05 800 

3PP0581C 

BPPC5B2G 

8PP05830 

3PPC5831 

BPP05832 

BPP0563S 

3PP05634 

3PP05F: 

3 PP C- . w 
8 PP •• 7 


SOURCE STATEMENT 


FORTRAN SOURCE LIST U8 3 


,./ . / 


X ( J » K + 9 )=; BSUKRN*AP( J,K)*TIME- ) 

X( J»K+ i 0 2 ) - A B $ f AKRW *AP ( J » K ) *T I ME ) 

GD TO r, '» 

v sh=v;jf ( j, K)/<( j,k + a6) 

SW= 1 . T-5Q 

V W i;: T = V W 

V N W = V 3 

S W L T = S W 
S N W = S 3 

I F ( G C • l T . G . 0 5 S W c T = SC 
I F ( GC. LT. ■. ) $NW=SW 
IF(GC.LT. , ) V N W = V W 
IF(GC. LI . . ) V WF T = VO 

rw5T*sxet*sw-;;t • 

RN W= Z • ) * VN W# S NW + SNW *S f i W* ( i . - 2. 0* V NW / V W ET } 

TGT= RWET/VWeT+RNrf/VNW 
RAW=Rd£T/VWc T/T3T 
RAN=RNW/VNW/TOT 
IF(SW.l-Q.SWET) RW = R AW 
IF(SW.EO.SWET) R3 = RAN 
1F(SW.£(J.SNW)RW*RAN 
IF(SW.EQ.SNW) RO=RA W 
X( J,K+94)=R0*V4 
X ( J*K+JLl»2)*RW*VA 
30 VOMAX = X ( J » K+48 ) 

vomi n= nn 

I F { GC • GT.O *9 ) V 0 MAX =,} . 96*X { J , K+48 ) 

IFCGC.LT }V3 MI N = u. 4*X( J,K+48) 
IFIVDMAX.LT.VQF ( J, K ) ) VQMAX=X ( J, K + 48 ) 

I F ( VOM I N . GT . VOF ( J , K ) } VOMI N=0. 0 
VWMAX=X( JtK+48) -VOMIN 
VWMIN=X(J*K+48) -V3MAX 
1F(FL3(J»K).GT.4*G)G0 TO 83 

IF(FL3(J,K).LU.4.2.GR.FL0(J,K}.EQ.3.2) GO TO 85 
V=VA+VQF ( J » K ) 

IF( V.3T.V0MAX) GO TO 82 
VOF(J, K)=VQF( J, K)+VA 


X( J,K + 1 IO)=0.C 
XI J,K+118)=VA 
IF(GC.GT.O„0)FLQ(J»K)=4.5 

GO TO ;,03 

32 CONTINUE 

X( + )=V-VCMA,< 

X ( J , K+ : 1 3 ) =VOMA X-VOF( J,K) 

; '.v'.tOFCJf^^VOMAX 

GO TO 'I" 10$ 

3 5 V=VA+X { J , K+43 ) - VDF ( J,K) 

I F ( V . GT* VWMAX ) GO TO 36 
V 0 F ( J » K ) = VGF C J t K 5 - V A 
X( J, K + i =V A 
X( J,K+il8 )=0.0 
GO TO 1 

36 CONTINUE 

X ( J , K+ 1 3 0 ) =VQF ( J t K ) -VOMI N 

, : . : .. ' :■ ' ; ■ ' 

' ■■■ ■■ ■ ■ ■., '^mSsM ■ sUS 8 



/XOOv; 

HaliifRiS® 


WtWMi 


siSS 


, C 

3 PPOE - .. -i , 
3 P P * <f, ; 

j-i p p 0 f 

6PP., 5 .*.:- 
3PP0tJb? 

bp Pi 5 c . 

dPP05d 70 
BPP05900 j 
3 P D 1 3 9 
3PP059C 
BPPO‘39 3 
8 P P 0 5 9 4 O \ 
3 P P C 6 9 
3PP0596- , 

i) P P Lc 9 7 
B P P G 5 9 o ; 
3PP0599 ! 

BPP06 ; 

B PP06C i 0 
BPP06020 ’ 

BPP06L -> 
BPP06-4 
BPP0605G 
BPP06070 | 

3PP06080 
BPPG6y9 ’ 

BP P 061 00 
3PP0611 0 
3PPC612-.; 
BPPQ6130 
8 PP06140 
3PPC&14i 
8PP0615 v* 

BP P 06 161) 

8PP06170 

3PPC61&4 

BPP0619 

BPPG6203 

BPP 06201 

3PPU>IA 

3 PP0622 J 

8 P P 06 2 3 j 

B PP062 4G 

3 PP 06241 

BPPO6250 

BPP0626G 

BPP06270 

BPP0628U 

BPPC629vj 

6PP063 0 T; ; 

3 PP 063 i 0 
3 PP 06 3 . ' | 

BPP ‘8; | 

i ? » G'.> v I 



ou»iv>-4-i>vjjps>#-*o-sionvji-P‘>— rv? /-*= 0^ \n -|> uj ro 


FORTRAN SOURCc LIST UB 

SOURCE S TAT ■: M F. <] T 

X ( JfK + i ‘.o ) = V4-X ( J,K+1, ) 

V 0 F ( J f K ) = V [ i M 1 N 
F L 0 ( J , K ) = 2 . 

GO T j 

3b V=VOF< J,K)-VOMIN 

IF(V.LT.X(J»K+94) ) G 0 TO £4 
VOF(J,K)=VOF( J»K)-X (J,K+9+) 

GO TO 9 * 

3 4 X ( J , K + 94 ) = VOF ( J , K ) - VOf ,1 N 
X{ J,K + 102 ) =VA-X { J, K+94) 

VOF ( J » K ) = VUMI N 
FLO( J,K)=2. t 
90 X ( J » K + A 1 0 } = X { J » K+94 ) 

X( J,K+L18)=X( J,K+102) 

P'J I F { C A S f * F 0 . 2 • i ) GO TO 15i 
X( J,93 )=Xl J,93) +X( J,K+94) 

X ( J » 94 ) = X { J » 94) +X ( J ,K + 10 2 ) 

GO TO 200 
150 CONTINUE 

X( J,K+94)=X( J.K+H ) 

X( J,K+i02)=X( J.K+lie! 

X( J ,83 )=X { J, 83) +X( J ,K+1U ) 

X{ J,84 )=X< J,84) +X{ J ,K + U8) 

IFUJ.LE.JR) X( J,K+24)=2.J 
IF(IJ.GT.I8) X{ J , K + 24 ) =2. 0 
I F ( X { J , K+ 24 ) .EQ.2.0) X ( J, K + 74 J =fl. 0 
IF(X(J,K + 24) .EQ.2.0 ) X{ J » K+84 ) = 3. 3 
I F I X { J,K + 24) .EQ.2.0 ) X t J , K+24 )=1 . 5 

3 X( J,K+il0)=0.0 

4 X ( J » K+- i 1 8 i 

5 230 CONTINUE 

7 RETURN 

3 END 

IBM4P ASSEMBLY UB 3 


ESSAG5S FOR ABOVE ASSEMBLY 


|0 7 


SOURCE STATcMlNT 


FORTRAN SOURCE LIST 



$1 BFTCSU H I 

SUbROliT I NE COND (KKK ) 

COMMDM/SYMARG/K ( 5 j * 16 >,AP(5G»8> , SC, VO, VW, TI ME , V3F ( 5 . t 
. XA ,FP(..' ,i! ), HYST?CAT, GAMMA, J,FLO(5G,S) ,GRID 
COMMON /CHIN/ 1C AiN{ 50 ,8 i,AWDIF(50,8),4DIF0(50,8),3AS.- 
GO TO (2? ,7 9), KKK 

C ENTER CHANGES IN OTHER CONNECTED NODE POINTS 

270 CONTINUE 
I) FORMAT ( 20X, 41 _> ) 

13 FORMAT (l NX » b Fi 3 . 5 ) 

N*IR*JR 
JIM = fN- JR 
VS=O.D 

1MBASE.LU.2. } VS = 1 .0 

DO 287 M = i,8 

1 F ( X ( J » M+ 2 4 ) • E Q • i • 5 ) GO TO 300 
IF( J1M.LT.1J ) 30 TO 273 
IF(X( J,M + 64) .LE.T. ! J GO TO 28 ) 

27 3 I F { X { J , M+ 2 4 ) • GE • 2 • ' ) GO TO 2 30 
300 CONTINUE 

I F ( X { J » M ) * r Q • j » >) GO TO 280 
NX=M+4 

IFINX.GT.8 ) NX® M-4 
MMK=I J+JR-M+4 
IF(KEQ.i) MMK=I J-JR+ 

IF(M.EQ*2)MMK=IJ+I 
IF(M.E0.6 ) MMK= I J-l 
IF ( M.EQ.7 ) MMK= I J-J R-l 
IF(M.EQ.B) MMK=I J-JR 
IF(MMK.GT.N) GO TO 28 
IF(MMK.LE.O) GO TO 280 
X ( MMK , NX + 56 ) =X{ J,M+56) 

XCMMK, NX + 1.6) =X( J » M + it) } 

IF(X( J,M+16) .£0.4.1 ) X< MMK,NX+i6)=4.1 
IF(X(J,M+i6) .60.4.2 ) X (MMK, NX+1.6 ) =4. 2 
IF(X(J,M+16).EQ.3.5) X ( MMK, NX+16 ) =3.0 
IF(Xi J,M+16) .EQ.3.0 I X(MMK, NX+16) =3.5 
I F ( X ( J i 73 ) .EQ.4.0) X( MMK , NX+24) =2 . 0 
IF{X{J,?3).EQ.4.0) X( J,M + 24)=2.0 
1F(X(MMK, 73) .60.4.0) GO TO 275 
I F < X ( J , M+24) .EQ.1.5 ) X( MMK , NX + 24) =2. 0 
IF(X(J,M+24). EQ.1.5) X< J,M+24) =2.0 
I F I X {MMK » NX + 24 ) .EQ. 2.0) X ( MMK, NX+74) =0. 0 
IF (X {MMK, NX + 24) .EQ. 2.0) X ( MMK , NX + 84) =0.5 
IF(X(MMK, NX+24) .FQ.2.P) X( J , M+74) =0.0 
IF (X{ MMK, NX+24) .FQ.2. 5 X ( J, M+84 ) =0. 3 
I F { X ( MMK * NX + 24) • EQ. 2. 0 ) GO TO 275 
274 IF{CAT.NE.,„. ) GO TO 275 

I F ( X ( J , M+24 ) . EQ. 1.5 >X ( MMK , NX + 24) =2 . C 
IF(X(J t M+24) .EQ.1.5 )X<J,M+24)=2.0 
IF(VS.NE.O.O) 30 TO 275 
IFCM.EQ.1.0R.M.GE.7) GO TO 272 
X ( MMK, 83 ) = X( MMK , 83 ) +X { J# M+ilO ) 

X (MMK,84)=X(MMK»84) +X( J»M+118) 

GO TO 275 


3 P P (j g A ?j 
3 P P • i 6 3 7 
) , I J , I R , JR , ti P P r Jn i b 

PPP06590 ; 

3 PPCoc. 

3 P P 5 6 _ 
BPP05O- 
8 PP 0663 j 
3 PP C66 4 . 
3PPC665 : 

3PP0665O 
B P P 0 6 3 7 J 
3PPG663 i 
BPP0669 
BP P 0670 • 

BPP067. (; 
3PP0673 
BPP0674', 
a PP 067 5.1 
3 PP 01 7 6 -■ 

3 PP 067 7 
BPP G67b 
BPP0679 • 
3PP068.. 
3PP0681 
B PPQ682 ' | 

B PP068iU . 
BPP0684 
3PP0685 . 
BPP0686^ 

2 PP Oo ' 7 
3PP0688 . 
3PP0689; 
BPP069: .. 
BPP0691G 
BPP0692 ■ 
BPP0693 
B PP06940 
BPP0695U 
3PP0696.J 
BPPC6970 
BPP 06980 
B PP06990 
BPP0700 J 
BPP,. 7 10 
BPP07020 
BPP07030 
BPP 07042 
BPP07050 
BPP07060 
BPP07=. 

B PP 07 * 

3PP0?:9 i 
BPP 07, | 

BPP.', 1 ! 


110 


SOURCE f 

> T A T 5 

272 CuNTINU' 


I F t X ( J , f 

i + 1 6 5 

IFtXtJ,? 

*1+ 1 / 

IFtXtJ, f 

4+1 is 


FORTRAN SOURCE LIST US 


:/ -/ 


• t Q • 3 • .0 A , X t J , M+ J 6 ) . E Q . 3- .5 ) GO TO 275 

) .ME. JC< PMK,NX+94) } AD I F3 ( J , M } =X { J, M+U }-X( 


, M K , T< + 9‘ 

(MmK, n x + 


IF(ADIF3(J,M).Mt . i .U.AND.AWDI F( J,M) .Nfc.C . 2 ) B ASe = i . J 
X ( MMK » NX + 9 h- ) ~Y. ( J.M+ll ) 

X ( MPK, MX+ I 2)=X{J,M+i 8) 

I F { BASE. FQ» 3 .0) I C A N ( J, M ) =MMK 

I F { Q A S c • i. U • , • ) P R I M T; G»X{J,M+74) , K ( J , M+ 64 ) ,XtJ,M+I6) , X ( MMK , NX + 7 


.XtMMK, NX+64) ,Xt MMK, NX + lb ) 

I F t BAS E . t-Q« 0 . C. AND • Xt J, M+24) 


I Ft BASE 


. « 0) PRI NT 


1MK, ICA 


) , X ( J » M + 94 5 » X ( J , M+ 102 ) 


0. AND.Xt J,M+24) .GE.2.0J PRINT 13 » X ( J ,M+i j 


M C J ■ M > 

) , X { J » M + 


CONTINUE 

CONTINUE 

GO TO 500 

FILL PROPER CONDUCTIVITIES 

00 78 K=i , 8 
B A S E = 0 . 0 

I F { X < J,K) .t g.0.0) GO TO 78 
IFtCAT.FQ. j . )*{ J,K+24)=1.0 

1 F { X < J,73) .tiQ.4.>*) X( J,73)=3.U 
SO = Xt J,K+56)/X( J , K+48 ) 

SW=1.0-S0 

IFtXtJ, K + 16) .fcQ.lO ) X( J,K+74)=X( J,K+32) 
IFtXt J, K + 16) .E0.2.G ) X < J , K + 74 ) =X ( J , K + 40 ) 
IFtXtJ, K + l 6) .EQ.l.O.OR.Xt J, K+16) .£Q. 2.0) 33 
IFtXtJ, K+16). NE. 4.5) GO TO 821 


-mm&m 


IFtXt J,K + 16) .E0.2.D ) X( J, K + 74) =X( J, K + 40) 

IFtXt J, K+16) .EQ.l.O.OR.Xt J, K+16). EQ.2.0) 33 TO 830 
IFtXtJ, K + 16). NE. 4.5) GO TO 821 
VL = VO 

IFt VW.GT.VL) VL = VW 

X( J,K+74)=X{ J, <+32 ) *VO/(V3*SO+VW*SW+6.0*VL*X(J,K) ) 

GO TO 830 

IFtXtJ, K + 16) .EQ.3.D .OR.Xt J ,K + 16 ) . EQ.3 .5 5 GO TO 322 

GO TO 823 

X{ J,K+74)=Xt J,<+32) *V0/( V3*SQ+VW*SW) 

GO TO 8?C» 

ARAT= VO/ V W 

IF (GC.GE.O.O ) GO TO 824 
IF(GC.LT.O.O) GO TO 825 
AKRW=$W*SW*X{ J, K+40 ) 

AKRNfcXt J,K+32)* (2.0*SO*ARAT+SO*SO*tl.O-2.O*ARAT) ) 
TOT = AKRW/VW+AKRN/VO 
X(J,K+74)=X{ J,<+3i)*VO*TOT 

GO TO B30 

AKRW^Xt j,K+-4 , )* (2.0*SW/AR4T+SO*SO*(1.0-2.0/ARAT) ) 
AKRN=K ( J,K+: '.'•:) *$0*$0 . 

T OT = A K RW / V W+ AKR H/VO ; 

X( J,K+74)=TQT*VD*X( J, K+3 2 ) 

CONTINUE 

X ( J , K + 74 } = AB S t X l J , K +74 ) ) ■ ; , 




X( J,K+i!O)»0.u 
X t J , K + 118 ) -§ .3 
X ( J . K +34 ) = . 

X ( J , 83 ) =0 . 0 


RPC'. Vv 
3PPC7 : 
BPP 07 t o 
CP Pi. 7 > 7 

3 P P 07 ... 0 
3 PP 07 .. 9 

, BPP :-i 

8PPC7 . 

5 PPQ'/i .. - 
33P W 07 
3 P P*,7 ? + 
BPP0725 
3 P P 0 7 2 6 J 
3 P P 0 7 2 7 
BP Pi jJl'i 
BPP0729 
Q P P 07 3 0 ' 

BPP i 
B PPG 7 0 2. 
BPP 07 3 - t. ■ 
B PP07340 
BPP 0735*. 

' BPP 07 So 
BPP 07370 
BPP 073 80 
BPPC739 
BPP 07 430 
8 PP0741 0 
3 PP074i .* 
BPPC743 > 
BPPC7443 
BPP 074 50 
3PP0746 
3PP0747 i 
BP POT 43? 
B PP 07490 
3PP0750C' 
8PP07510 
. 3PP07520 
BPP07530 
BPP4S750i 
BPP 07347 
BPP07550 
BPP07560 
' 3PP0757C 
BPP07571 
BPP 07 580 
B PP07 590 
BPPC76 ) t 
BPP076’ 

■ 3 P P C ? 6 ; ...> 

JP’PC I 




ili 


FORTRAN SOURCE LIST U6 4 

‘ / -t 

>N 

SOURCE STATEMENT 


:1 

X{ J,64)=0 . j 

3 P P 0 > “ ' 


IK=J/JR+1 

BPP 07 65 

.3 

JK=J.-JR*{ IK-U 

h a i .. .. : 

:4 

IF ( JK. LE . i ) 1 K= 1 K _ i 

3 PP076 1 u 

'7 

IF(JK.LC.O) JK=JR 

3 P P C 7 * . ; 

#L 

IF ( IK. i:Q. i )GO TO 2.4 0 

3 PP07 69 

5 

IF( IK.rU.2.AUL>. K. ;:Q .1 )G0 TO 140 

& R D 7 7 1 


I F ( IK. c 3 . ... . AND. K . GE .7 ) GO TO i 40 

3PPQ77 0 

■3 

R=GC* ( ( » 0 / X ft ) — ( ; . rj/Kl J,K) ) ) 

3 PPG 7/.. 

■4 

R = - R 

3 p p . j rrs 

■5 

J U T J ' 

3 P P 0 7 7 t vi ; 

■6 1 40 

R = G C/K { J, K) 

3PP0775C ! 

•7 150 

CONTINUE 

3 P P L 7 7 C 1 

0 

IF(X{JjK+a.6) .EQ.3. ) X ( J » K + 8 4 ) = R 

B P D C 7 7 7 

3 

I F ( X ( J » K + 16 ) . EG. 3. 5 ) X ( J j K+84 ) = -R 

6PP07 7 . 

6 78 

CONTINUE 

3 PP 07 V 9 vi 

D 500 

RETURN 

3PP<J7n 1 

1 

LND 

B P P 07 bi, . : 


I3MAP ASSEMBLY UB 4 

i Z / 0 ■; L 

ESS4GHS 

FOR ABOVE AS Sl.MBLY 

'■ S 
1 


FORTRAN SGURI 


LIST 


Xi 


sour; 


STATE Mt NT 


$1 BFTCSUB 5 

SUBROUTINE CHANGE 

COMMDN/SYM ARG/ X ( 5 Ot 160) ,4P( 50, 8) , 

IX A ,FP{ . J, 11), HYST, CAT, GAMMA, J, FL3 1 5C , 8 ) . b 


6 9 


6 80 


C 

c 


55 


u 

3 

6 

7 

0 

1 

4 
? 
2 

5 

6 
7 

1 

2 

5 

6 
7 
2 
A 

5 

6 

7 

2 

3 

:4 

|5 

b 


55 5 
556 


: f » . V ■# j ' * ■ i f > « » * 

COMMON/CHAfi/DIFD,WD IF 

" " ' , 10), I, I CHANG, AR AT, TRAP 

\ t *. a t r* n 


3C , VO , VW, TI ME , VQF { 50 , 8 3 , 
RID 


I J , I R t J 


, I CP 


COMMON/SOLB/Cm? ( . 0 
CUMM'j 4/SOL J iOO,CP, JIM 
DIMENSION 3 ( ) , C ( b ) 

R A L K a 

, 5Hu! F3=, El: .5,5X, 5HWD I F = , £13. 5 ) 
AND. I .GE. .7. AND.J.EQ.22) PRINT68 
AND. 1 .GL. . 7.AND.J.EQ.22)PRINT63 
'ID. J.EQ. 22) PRINT5 3 


A ■ 


FORMAT ( 

1 H I CP. 

I F( I CP. 
i f u : ? . 

I F ( ICP.FQ 
FORMAT { IX, 


IQ, 

; 'U ■ 
*Q ■ 


. A N 


,01 FO , WD I F 

, (X( J,K) ,K = 17, 
(XIJ,K) ,K=?5, 


DO 

Ur 


553 


2. AND. I . GE... 7 . . . . _ 

2.AND.I.GE. ! 7 . A N D « J .FQ.22)PRI MT680, ( X { J , K+94 ) , K= 1 
BE *5.7) 

,8 

FOR FLOW CHANGES SO 
IN EACH NODE OF THE 


554 


555 


K = 

select tubes 

IS SATISFIED 
C { K ) = D . 0 

CON 1INUE 

ETTIMS FLUID IMBIBITION 
BASE=0.0 

« a =c.:. 

CAT= ,. ■ " ; ■' 

AC F — 0 . ! ) ■ 

IF(GC.GE.O.D.AND.WDIF.GT.O.O) 
I F ( GC. L T . ' . ! * . AN D . D I FO.GT.O.O) 
I F ( K A . fc Q . >> . 0 ) K A = 4. 5 
IF ( KA. EQ. 4.5 ) SAT=0.3 
B ( 1 ) = D, . 0 
DO 555 K= 1 , 6 

v=o.o 

I F ( X { J , K + 2 4 ) .GE 
IF(C{K).NE.O.G) 

I F ( FLO { J,K).Nfc 
I F ( AP ( J , K 5 .LE. 

V = AP( J ,K )*;■',( J, <+7 

V = A B S ( V ) 

b ( K P { K 3 + V 

CONTINUE 
IF(B(9) , t;Q. 0.0) 

RN=RNDY1 U> . ) 

DO 557 K = E , 9 

\ ' I F ( RN • L 7 . B ( K ) / 3 ( 9 ) ) GOTO 558 

65 T-. CONTINUE 
55 8 K«'K— 1 

GO TO '.570 

PLHDJLAR FLOW CHANGE 
I F ( CAT. GE . . J) GO- TO 560 
KA=4.5 

CAT = . • . 3 
GO TO 554 

I F { C A T . E Q . ; . . 3 ) GO TO 
IF(CAT.EQ. N43 GO TO 


THAT MATERIAL 
MODEL 


BALANCE FOR EACH 


KA=1 .0 
K A = 2 . 0 


.2.0 )C 
GO 13 
GO 
GO 


. KA ) 


TO 

355 

TO 555 
TO 555 


) 


GO TO 559 


55 9 


55 


561 

562 


BPP '!?■: ■: 

3 PPC7 . <■ 

. , 3 P P C7 . + : 
3 PPG 7 3 5 
ti P P 0 7 6 c> 
3PP0V 3 5 1 
8 P P t» 7 3 £ g 
3 P P C 7 87 
BP PC 7 , i 
3 P P 07' 3 9 j 
3 PP 07 3 9,; 

B PPG 7 3 9 4 
BPPC7 69. 

3 PP'07 8 94 
3 PP €7 695 
BPP079 
PHBPP079 ,,u 
BPPOlyCC 
BPP079-. 

BP P 07 9 4 
8 P POT 9 50 
8PP0796U 
3 PP 0797 
BPP0793 
B PP 07 98 i 
3 PP 07 990 
3PPC8 . 
BPPC8U1 
3PPOBOZO 
BPPCbwi • 
BPP 0894,4 
BPP0835' 1 
3 PP08060 
3PP08 V 7 
BPPQ8..8 
BPPGd/Bi 
BPP08090 
3PP081 j, 
BPPGBii 
3PP08123 
3PP0813U 
BPPQ314: 
6PPC815C 
BPP031 60 
3 PP 0817 0 
3 PP 08 130 
BPP0319C 
BPP08200 
3 PP08210 
3 PP 08 223 
BP P 0825 . 
B PP Ob 2^-0 
3 P P 08 ? 5 
BP! * ! 

B ? P -■ ; : : 7 



113 


SOURCE STATE Mt NT 


FORTRAN SOURCE LIST UB 5 



I F i CAT . EQ . .5) GO TO 563 
IF(CAT . L'Q. C . 6 ) GO TO 564 
I F ( CAT . EQ. G. 8 >30 TO 565 
GO TO 620 

56 1 CONTI NUF 

C SLUG FLOW CHANGE 

I F ( WDI F . GT . u . 0 ) K A =4 .2 
IFtDIFO.GT. K A = 4 • 2 

C A T = 

GO TJ 6 54 

5 61 1 F ( W 01 F • 37 . 0 • 0 ) K A = 4 • l 

IF(DIF;j.GT.0.,JKA = 4.i 
C A T = ,, . 2 
GO TO 5 34 

562 CONTINUE 

C DISPLACEMENT FLOW CHANGE 

IFCDIFU.GT.G. J) KA= 3. 

IF{ WDI F. GT .0.0) K A = 3 . 3 
CAT*!). 5 
GO TO 554 

563 CONTINUE 

C NON WETTING FLUID IMBIBITION 

K A = Q . D 

I F ( C I F LI • GT • C . 0 ) K A = 2 . 1, 

I F ( WD I F . GT . ,j • 0 ) KA = i 
CAT = : .6 
GO TJ 554 

564 CONTINUE 

I F t WDI F . GT . 0 . 0 ) K A=3 .0 
IF{ DIFQ.GT.0.j)KA=3.5 
CAT=i). 7 

IFIKA.fcQ. 0.0)00 TO 620 
GO TO 554 

57o V=AP(JtK)*X( J » < + ?**) *T I Mg 
FLOU, K) = KA 

IF (KA. 00. 3.0) FL 0 ( J» K)*3 .0 
IFIKA. ,0.3.5 )FLO( J, K)=3.5 
I F { K A . L U . 4 . ; ) FL 0 ( J , K 5 =4 . 1 
IFIKA. E Q.4. 2)FL0(Ji K ) =4 . 2 
I F { K A . £ 3 . * . 0 ) FLO ( J ,K )=_,.Q 
IF(KA. _Q.2.u> FL0( J ,K. )=3.5 
IFCCAT.hU. .5JFL01 J,K)=4.5 
VOMAX=X( J,K+48) 

V0MIN=0. 0 

"I F(GC. LT. . ) y/OM I N =, . 04*VQMA X 
’ IFCGC.GE.O. ) V0MAX=2.96*V0MAX 
IF(Xt JfK+56 ) . GT . VO MAX } VOMAX=X( J,K+48) 
IFIXI J,K+56) .LT WOMIN) V3MIN=G..O 
VWMAX=XU»K+48)-V0MIN 
VWM IN* X ( J » K+48 ) -V3M AX 
VOLO = VOMAX-X ( J* K + 56 ) 

VOLW=X ( J, K+56J-V0MI N 
CA SE = 2 * 0 

IFiDIFO.GT .1.0) V=X ( J » K+102 ) 
IFUDIF.GT.D.O) V=XU*K+94> 


3 p D (j C r, 

3 jpo - • . : 

3 P P i F . , | 

BPPOF \ 
UP POP . < 

B PPG.' 3 .. 6 
•3 PP LB '■ - : 

BPPLF-:.i ; 
BPPOo. 4 0 . 
BP°C655 j | 

3 PP 06 0 f. » 
BPPut^G-: 1 
BPP 03 . 6 
3 P P 083 c4 ; 
B P P U ' 1 S o 5 
BPPC837. i 
BPPO'.'.GO ! 
3PP08390 
3 PP Cot 
B P P 08 4* .■ | 
3 P P O’’. *. 2 u | 
BPP0S4-. i 
B PP0&44, f 
3PPCE45 
BPP0846) 
BPPG847 , 
BPPG849.; 

B PP 08491 . 
3PP0P402 i 
3PPG8493 
BPPG8494 
3 PP 08495 
BPPG8496 
BPPCS5 
BPPD65) „ j 
BPPG852U , 
3PPC853 _ 
BPP0854.r 
BPP 03550 
BPP08560 
3PPC35 7i/ 
BPP 08571 
B PPQ8 580 
BPP06590 
3PP0B60U 
BPP0361 j 
SPP0862G 
BPP 08630 
3 PP0864U 
BPP 08 650 
8 PP0B66 0 
3PP08670 ^ 
BPP 08 6 c . | 
RP n <' : ;-jG I 
B PPG c TOO i 



li 4 


SOURCE STATEMENT 


FORTRAN SOURCE LIST 118 5 



IF( WDIF.GT.G.O) G3 TO 601 
I F ( DIFO.GT.U.G) 30 TO 603 
C WATER EXCLSS 

601 IF(V.GT.WDIF) 3D TO 6 2 
WD I F=W ul F-V 
X( J» <+9 4) = j.v 
X( J,K+; C)=v! j,K+i 2)+V 

IF(V.LT.VOLW) GO TO X. 

FLO ( J» K )=2. j 
VOF (Jr K. ) = VGM J N 
GO TO 35% 

5 . VOF(J}K!=X(JiK+36)-V 

GO TO X 2 

632 VOF ( Jt K)=VOMIN 

I F ( V3LW.GT.WDIF) Vu F { J, K ) =X ( J » K + 56 ) -WDI F 
IFIVrjLW.LT.WOIF) FLO( J , K ) =2. 3 

IF(GC.LT.O.O.AND.FLO( J,K> .EQ.2.0) FLQ( J,K)=4.5 
X( J, K+94) =X( J » < +94) -WOIF 
X ( J »K + 1 1 2 ) =X { J i K + l } 2) +WD I F 
WDI F = 3 

12 CONTINUE 
SO=V3F(J,K)/X< J,K+48) 

IF(FL3( J,K).£Q. 4. 5 .AND.GC.GT .0.0 ) CAS5=C.5 
IF ( CAS E.EQ.O. 5. AND.SQ.LT. 0.80) FLO! J,K )=4.1 
1 F C FLQ( J,K).tQ.4.2. CR.FLOI J*K).EQ.4.1 ) BASE=i.O 
IFtBASe.tU.l.C. AND. GC.GE.i .0) CASE =1.3 
I F (BAS E. EQ.l.L. AND.GC.LT. u.O) CASE =2.0 
1F(CASE.EQ.*.0. AND. SO. G£. 0.8) FLO( J,K )=4.5 
IF( CASE. EQ.2.0. AND. S3.LE.0.2 ) FLO ( J , K ) =4. 5 
CASE=3.0 
BASE=3.:> 

IF( WDIF.NE.O.O) GO TO 3 54 
GO T3 590 
C OIL EXCSSS 

633 I F { V.GT.OIFO) 30 T3 6 4 
DIF 0=0 1 FCS-V 

X( J,K + 94)=X( Jr < + 9 c i ) +V 
X ( J » K + i. . 2 ) = • . 

IF < V.LT.VOLO) GO TO 1., 

VOF { J j K ) = VOM AX 
FLO ( Jr K)=I, . 

I F ( GC . GT . . ) FL0(J»K)=4.5 

GO TO £ 5 4 

13 VOF(J,K)=X( J,K+56)+V 

• - 1 fin th • c, 

654 VOF ( J» K ) = VUM AX 

IFIVPLQ.GT.D1F3 ) VOF( J , K ) =X ( J , K + 5 6 ) + DI FO 
IFtV0L0.LT.DlF3 ) FLC( JrK) =1.0 

I F ( GC. GE.O .1 .AND. FLO( Jr K) .EQ.l.U) FLO t J,K)=4.5 
X( J t K+102)=X(J,K+1 72)-DIFQ 
X( J,K+94)=X( J,K+94) +D1F0 
DIF0=3.0 
4 CONTINUE 

S0=V0F(JjK)/X(J»K+48) 

I F { FL3 ( J,K ).F:Q.4.5. ANO.GC. LE.0.0 )CASE = 0. 5 


3PP0CT.t 3 
BPP'.E/'.. , 

3 p o /3 

yp p uS7 ^ 
dPPQt i 5 . 

J <-> P r p. 1 1 
B PP C *. i 7 
B P P *7 & 1 < > ; 

d P P 0 \ / 9 0 

bp •, 

BPPC93. 

J P P Oo o j ! 
B PP08F40 . 

3 ppgbpS : 
B P P 8 o 
BPPOBb fj > 
3PP0 ! . ,4 

B P P u 3 8 9 1 > 
BPPOBQ 
a p P 0’ : o ; 

3 PP Ob 92 0 J 
3 PROS 93 
BMP OR 94 j 
3 PP 08950 ■ 
BPPQ&96.) 
BPP0897' 
BPPC398 » 

8 PP08990 . 
3PPt 9 
3PP09 *1, 
BPP0932" 

3 PP 0903 0 
3PPU9 + 
3PPU9 <5 
8PP 09060 
3 PP 09070 
3PP09- J ; 

BP P 09 9 0 
B PP09100 
B PP 091 i 0 
BPP0912 / 
BPPuVj 3 # 
8PP09: 40 
3 PP 091 50 
3PP09i6 
B P P 09 1 7 4 
BPP09130 
3PPC9190 
BP P €9203 
BPP092I j 
3 PP 092 2 0 
8PPC923 
b P P09 ** 
BPPf 9 , 3 
p H b 09 . 



SOURCE STATEMENT 


FORTRAN SOURCE LIST UB 5 


.. Z / 


IF (CASE. CO. 0.5. 4N0. SO.GT.3.2) FLO( J ,K )=4.2 
IFIFL3(J,K).EQ.4.2. OR . FL3 { J , K ) . SQ. 4. I ) 8 ASE = _ . 
IFIBASF.EO.i.C. AND. GC.LT.U.O) CASE *1.3 
IFIBASE.FQ.l.D. AMD. GC .GE.J.i)) CASE-2. 3 
IFfCASF.EQ. ' .0. A,ND. S0.LE.D.2! FLO{ J,K)=4.5 
IFICaSF.LQ.E.v.hND. SD.G&.0.8) FLO{ J » < ) = 4. 5 
BA SE= , . 

• A o Z — * 

1 F ( D 1 F j . N S » r J . 0 ) 33 TO 554 

i9o conn no-:.: 

6? IFCDIFT.Ll . • * * i ; 1 )DxF 3= .' 

IFlWDIF.Lb. .j'fc-1 ) W 0 1 F = 3 . 

1F{ DIF 3.: 0.0.0. AND. WOI F. PU.Q . 0) GO TO 633 
ACc = AC + 1 . 

DO 35 K=i»S 
35 C(K)=D , 

CAT =0. i 
K. A=4. 5 

I F ( A C E . E Q . 1 . ) GO TO 554 

63 ' IF( WOIF.NB. 3.0. OR. DIFJ.ME.0.0 ) PRINTS 90, DI F3.HDIF 
RETURN 

END 

FORTRAN PROGRAM UB 5 

t FUNCTION SUBPROGRAM REFERENCES 

! ' ' • ■ • ' . 

I : 

IBMAP ASSEMBLY UB 5 


sppon;. . 

3 P P C 9 !i 3 
6 PP 99 2.9 

3 P P 39 i 0 j j 
3PPQ9j i 0 1 
3 P P ■ £ . . 

6PPQ933 

B ? P 09 3 ; 

6 PP 09.-5 ju 
3 P P 09 s 6 

g p p {jq 7 . 

BPP09;:,.o 

B HP 093 px 
3PP093S: 1 
BPPC9 jop ■ 
a P P 09 > J H 
3 P P 0 9 P 3 , 
3 PP 09 386 
B P P 09 3 3 7 
BPP09390 ! 
3 PP 094 j „ f 
BPP0941 '■ \ 
12/.' 

1 


1 If - 


IESSAGES FOR ABOVE ASSEMBLY 



SOURCE STATEMENT 


FORTRAN SOURCE LIST 


0/ -I 


>1 3F TCSUB 6 

SUBROUTINE PDIF 

COMM3N/SYMARG/X (52, 16 ),AP(59,8) , GC» VO, VW , T I HE , VOF t 5 
,FPU0,10) T HY ST, CAT, GAMMA, J, FLO (50,8) , GRID 
PRESS DIFFERENTIALS FOR ALL TUBES AT EACH MODE POINT 


» 8 ) , i J 


1.XA 

FILL 
J = IJ 
IK=IJ/JR+ 
JK= 1J-JR 


1 F ( J K . L 


UK- ) 

) IK =I<-. 

I F( JK.Lr. ) J K = J R 
F0 P M .,T ( i X * . 1 ) 

M J,?4) = l.P (IK, JK) 

I F ( I K. . Q • ’ ) G 0 TO 3 7 

1F( JK.lQ.I JGU TO 9 c 
I Ft IK.,:Q.1R) 03 T"J 96 
IF(JK.EQ.JR) GO TO 10u 
X { J , 65 ) =EP 1 1 K, J K } -'£ P( IK 
X t J,66 ) = t:P( IK, JK)-PP( IK, JK 
X(J f 67 ) =E P ( I K , J K ) ~fc P U K+ 1 , 
X( J,63)«EP(IK, JO-:P(IK+l> 


i, JK+1) 
JK + 1 ) 
JK+I ) 


3 P P 6 9 -i- . 
hpP-.OA. 

I R , J 3 , 6 P P CO -r '•» J 
B P P 094 0 ■.< 
3 P P 0 9 6 
6 P P U9 4 7 

JPPO >• 

3 PP 09 A 00 
3PP KO 
BPP UU 
3 P P 09 5 1. j 
■3 p P 0 9 3 
3 P P 0’° 5 4 
P P J9 5 5 
8 P P 09 56 U 
5 HP 09. 7 
3PP.0US 
BPP0959 
B PP 09600 
3PPU9U. 


X( J,69)=£P(IK, JK > — E PC IK+1, JK-1 ) 
Xt J,73 )=EP(IK, JK)-FPt IK, JK-1 ) 
X(J»71)=£P(IK,JK)-E P( IK-1, JK-S ) 
Xt J,72)=EP( IK, JK )-EP( IK-1, JK) 

GO TO 132 

37 IF( JK.ea.l ) GO TO 88 
IFtJK.EQ.JR) GO TO 90 
X t J , 65 ) =0 • 0 
X ( J , 66 5 = 0 . 

X { J,73 )=(>. 

Xt J,7i )=0.0 
X t J , 7 2 ) =■ • . ■ 

Xt J,67)=EP(IK, JK)-EP( iK+1, JK+4 ) 
X{J,68)=EP(IK,JK)- P ( i K + * » J K ) 

Xt J,69)=EPt IK, JK )-EP( iK+1, JK-1) 

GO TO rz 

3 8 Xt J,67)=UM IK, JK)-': PI kK+l , JK + 1) 
DO 89 K=l ,C 
IF(K.lQ.3) GO TO 39 
X ( J , K + 6 4 ) = . 

39 CONTINUE 

GO TO A >2 

90 Xt J,69)=tP(IK, JK)-EP( IK+1, JK-1) 

DO 9., K = ... ,0 

V'v 1 ;/ IFtK.EQ.5) GO TO 91 
X t J » K+ 6 4 ) = . 

91 CONTINUE 
GO TO 132 

9:. IFtIK.E3.IR) GO TO 94 

X( J,65)=£P(IK, JK)-FP(1K-1, JK+1) 
Xt J ,66 )=EP( IK, JK >-EP(IK, JK+1 ) 

X t J , 67 ) =E P ( I K, J K ) -E P( IK+1 , JK + 1 ) 
DO 93 K = 4 , 8 
Xt J»K+64)*0.Q 
93 CONTINUE 


3PPC963 I 
5 P P 09 6*+ J 
3PPC965 * 
BPP0966 > 
BPP 09 670 
BPP09680 
3 PP09693 
BPPQ97 > > 
BPP0971 0 
3PP09720 
BPP0973E. 
BPP09740 
BPP0976 . 
BPP0976C 
BPP0977 } 
BPP09783 
BPP09790 
BPP09800 
BPP C981y 
BPP0982 5 
BPP09820 I 
3PP09840 
BPP 09850 
BPPG986C 
BPP09870 
3 PP 09880 
BPP09890 
BPPC993J 
3PP099i0 
3 PP C992 ■; 
8PP0993E 
8PP099VJ 
BPPG99C u 
.DP '9 96 



1/7 




FORTRAN 

SOURCE LIST UB S 

.1/ 



SOURCE STATEMENT 





GO TO 132 


3PP 

7 

94 

X{ J,65 )=f'P!IK, JK > - T P ( iK-l, JK+1J 


b i-‘ P *. '} i 

2 


DO 95 K=2»8 


BPPCN993 

2! 


X( J,K+64)=0.C 


BPP10 

22 

95 

CONTINUE 


3 P P J ' 

?4 


GO T3 1 2 


PPP1 2 

15 

9 6 

IF( JK.CS.JR) GO TO 96 


B P P if. 

57 


X ( J » 5 5 ) = = P { I K , J K ) - : P U K-l , JK+ 1 ) 


3 P P 1 •. »' , > 

11 


>: ( J , 77. ) = tP ( I K, J K ) -t P{ IK-l , JK- ,1 ) 


3 P P 1 5 . 



X ( J * 7 . ) = ; P ( i K , J K ) - : P U K- i f J K ) 


3 P P 1 6 

B 


00 9? K ~2. t b 


bp p j.' 

B 


K( J,K+{>' ) = 0 • 3 


3 P P 1 0 J h u 

55 

97 

CONTI NUc 


3 P P i 9 

57 


!-■? ( J T [J x 2 


' K P P 1 , 

, 

9 8 

XC J,7* )=tP(IK,JK)-EP( IK-l, JK-i) 


B P P i [< i ; '■< 

•i 


DO 99 K = 1 , c 


3 P P 1 2 



1 F { K . E 0 . 7 ) GO TO 99 


3PP 1' 



X ( J , K+ 64 ) = . 


BPPI 7 4 - 


99 

CONTINUE 


b PPlOx 5.j 

50 


GO TO 102 


3 PP 1 . 6 


l J V. 

X ( J , 69 ) =EP ( I K, JK) -F PU K + l , JK-1 ) 


3 P P 1 j ( 

52 


X{ J,7*> ) = EP(IK, JK)-f-P{IK,JK-l) 


BPPI -xi, . 

; 3 


X( J,71 )=EP(IK, JK)-E P( iK-l, JK-1) 


BPP10190 

54 


X( J»72)=0.0 


3PP1G2-/..; 

I 5 


00 m K= 1 j 4 


BPPI *21 

56 


X( J,K+64)=< . 


BPP10220 

57 

131 

continue: 


BPP10230 

51 

132 

CONTINUE 


3PP 1024s.* 

C 


FILL POSITIVE DIFFERENTIAL PR. AND 

CAP. PR. 

RPP 1"'25. 

52 


DO 1.4 K=I *6 


BPP10260 

,3 


VOF ( J » K ) =0 .0 


BPP10270 



API J, 0=3. 


BPP10283 



I F ( X { J » K ) » . .. . ) GC TO i : 4 


BPP13290 



X( J,K + 64 )=X{ J f K + 64) +X ( J , K + SM 


BPP10300 



IF(X(J,K + 64) .GT.0.0 )AP( J,K) = X( J,K + 64 ) 

3 PP 1031 0 


134 

CONTINUE 


B PP 1X323 



RETURN 


B PP 10333 



END 


B PP 10340 



IBMAP ASSEMBLY UB 6 

; A / 

jESSAGES 

FOR ABOVE ASSEMBLY 




IIS 


FORTRAN SOURCE LIST 1 / 4 

SOURCE STATEMENT 


$1 EF TC SUB 7 

3 P P 1 0 i 3 'j 


SUBROUTINE SIMU 

3 P P 1 3 i 

C 

PROGRAM SIMU 

BPPi, i V 


COMMON/ SYMARG/X { 51 , 16 ), API 50, 8) , GC, VO, YW, TIME, VjF ( 5 ,8) 

, IJ,IK, JRfBPPl 


I X A , E p r , ,1 3 , HY ST, CAT, GAMMA, J,FLO{ 50,81 , GRID 

BPP 1C 90 


COMMON/ SOLC/AuQ, CP, JIM, ICP 

3 PP IDt j > 


COMMON / S I M/COM ( .5) , RA U ( j.5 V » KL *LK 

3 P P 1 v. 


DIMENSION F( 5 , 3) 

BPPi 


FACT=S ■ S A 7 . * .15/., .7 32 

3 P?lu- ■ u 


DC 5.' J= : , JR 

3 P P i f “ -* 


DO 6 I ~ , i R 

BPPI 4 5 


1 S K = 1 1 - , 3 * J R + J 

3 P P 1 . <rb - 


DO 20 m-i, e 

HP Pin A 7 C 


X { I S K , M + 8 ) =1.0 

3 PP 1 ,4b:. 


X { I SK, M + l 6 3 = .5 

3 HP! 49 


X ( I S K , M + 2 A ) = 1 • *5 

3 P P 1 3 


X { I S K » M + 3 2 3 = , : . • D 

■ BPP 1C 5., 0 


X ( I SK, M + 40 ) = ) .0 

8 HP I .-5 2 


xi isk, m + '+ 8 ) = ; . 5 

BPPI 53 


X ( ISK, M + 5 6 ) = . 3 

BPPlDS-u 


X( ISK, M + 6* ) = i . 3 

BPPIOOSU 

20 

CONTINUE 

BPPI'. 56 ■ 


X( ISK, 733=1. 

BPPi:) 57 j 


X{ ISK, 74)=a. 

BPP10580 


IF(J.EQ.JR) GO TO 1 20 

B PP10590 


DO 70 K = 1 , 4 

3 PP 10 6 C’i/ 


RN= RND YJ. ( Y 3 

BPP 1261 'J 


DO 75 M-l , LK 

BPP1D620 


IFtRN.LE.C DM (M3) GO Ti *08 

3PP1063C 

75 

CONTINUE 

BPP 10642 

ill 

X ( 1 SK, K ) =R AD ( M ) 

BPP 10653 

T 

CONTINUE 

SPP10660 


GO TO 130 

SPP10670 

1 2 0 

DO 1 60 K — * , 4 

BPP 1968%* 


X ( I S K , 8 ) = < . 

BPP 11/69 2 

160 

CONTINUE 

B PP 10700 


X (ISK, 3=0. ... 

3 PP 10712 

1 3 ' 

IF(J.N-.i) GO TO 5 v 

BPP 11)7 20 


DO 1 65 K=-,C 

BPP 10730 


X( ISK, K 3 =;!.-• 

B PP 10740 

155 

CONTI Nil . 

3 PP 1 75 

3,5'., 

1 F ( I . N . I 3 GO TO 1.5 5 

BPP 1076 u 


XIISK, ..)= . 

BPP10770 


X( ISK, 2 3=0.0 

BPP10780 

■ :; T: 

X( ISK, 6 3 =0.0 

BPP 1079 0 


X ( I SK, 7 3 => . 

BPP 1080 i 


X ( I SK, 8 3= . 

BPPU'taO 

155 

IFU.NE.IR) GO TO 180 

BPP 10 820 


DO 195 K=2 , 6 

BPP 15 83') 


X( ISK, K)»0.y 

BPP 17 34 J 

19 5 

CONTINUE *•«-, 

BPP 10850 

C 

FILL ALL TUBE RADII 

BPP 10 860 

1 3 0 

IF(J.EQ.i) GO TO 190 

BPPI 87 . 


DO 135 N=i ,3 

BPPI 9 


IF( I.EQ.l.ANO.N.NE.l) GO TO 115 

BPP 10 ,91 


i 


Ilf 


r 

'■w> 




,f ; FORTRAN SOURCE LIST Ub 7 


/ 


IF( I.EO.IR.AND.N.NE .3.1 GD TO 125 
I P=I -N+2 

I $=(IP-; ) * JR+J- ] 

X ( I S K , N + 4 ) = X { I S , N ) 

135 CONTINUE 

IF(I.tQ.I) GO T O j 9 C 
1 S= U-2 )* J n + J 
XtlSK, ,5 } = X ( I S » 4 ) 

*9U CONTINUE 

C FILL LENGTH IN 9-15 

DO 2S ) K= 9 , i, h 

XI 1 S '< 5 9 ) = '.,47.42 

IF {.XU SK,K-. ).=t).i .0) X(I3K»K)=1.4142 

230 CONTI NHL 

C FILL" OIL AND WA TER CONDUCTIVITIES IN 33-406 ND4* -46RL SP F.CT1 V^L Y 

00 3 0 . K = i sB 

X( ISK.K + 8 >=X( ISK,K+8)*GRID 
6 ( I SK, K)= { FACT* < ( ISK, K) **v) / X ( ISK, K+8 ) 

X( ISK, K + 3 2 ) = R ( I SK » K )/V0 
X( ISK, K + 4, 5 =B ( I SK , K )/VW 
C COMPUTE VOLUMES OF TUNES 

X ( I SK, K + 48)=3.‘)*1 ,73?*X( 1 SK,K+8)*X ( ISK,K)**2. 

C FILL FLOW OF OIL IN ALL TUBES 75-64 

X{ ISK,K + 56 )=XUSK,K+4U 
X(ISK»K+74)=X{ ISK, K +3 2) 

X ( I SK, K+84 ) =' .0 
X { I SK, K+94 ) = '.T 
X { I SK, K + lOL ) =0. 0 
X( ISK,K+1*0)=U.0 
X ( I SK , K + 1 i 8 ) —U « ./ 

X (ISK, K+128) =0. > 

X{ ISK, K + ? 36}= . . 

X { ISK, K+144) = :. 

X( ISK, K+ 152) =0.0 

IF(X(ISK,K). Q. . ) X ( I S K » K + 2 4 ) = 2 . 0 

333 CONTINUE 

X C I SK, gr )= . V; * . 

X{ ISK, 64)- 0 . 0 . 

X ( I SK , 9? ) = . 

X ( I SK, '-J '■■■) = . 

X( ISK, 127)-. 

X( ISK, I 29) =0,0 

C ASSIGN VALUES TO DEL PR IN 65-72 

IAN* IR-i 

'N;U* AN=IAN 
'• ! ' 1 !f?l^*CP/AN 

CO ‘301 K=3 , 5 

X C ISK, K+64 ) = Z 

Uii. CONTINUE 

X(ISK»71)=-Z 
X( ISK, 72 ) =-Z 
X { I SK, 65 ) =-Z 
X{ ISK, 66) =0.0 
X( ISK, 70) =0.0 
C ASSIGN PRESSURE 74 


OH P 10901 
3 pp i : 9 : ; 

B P 0 1 9 
l>pp 1 93 

p. P P 1 o 9 -v ; i 

3 P P i ' 93 
3 P ? 1 9 a 
BPP1 9 7 ; 

PPPIuO.U' 

3 P P JO 9 u ■ 
f. P ° 1. 
dHP I* 4/ 
BPP 11 u.' G 1 
3FPII - 3 
3PP 1* . 4 
B P P 1 1 50 I 

3 PP1, 06 C I 
3 PP 1" ,?o I 
Bppii s ; 
BPP ll 0 9u ! 
3 PP ll.i 00 I 
BPPII ix 
BPPlli2> 
BPP1M30 
3PP13 140 
BPP1U5 ; 
BPP1U6,; | 
BPPII. 170 
3 P P 1 1 B ’ 
BPP 11x9;, 

B PP 11 20 - 
BPP 1) 20 1 • 
BPP 112 2 0 
BPP 11 230 
3 P P 11 24 ;■ 

BPP lx 23 i 
BPP 11260 
BPP 11 2 70 
8 PP 1126 
BPPII 29 * 
BPPII 300 
3 PP 11 3 10 
3PP113z„t 
6 P P 1 1 3 3 * 

0 PP 11 3 40 
3 PP 11.35} 
BPP 11360 
B PP'll 37u 
BPP 1.13 3 0 
BPP li 39 . 
BPP 114 ; 
BPPuX 0 
B PP 1,1 4 ; 




5 p p 1 
3P.P1 
3 P P i 
H P P 1 
3PP I 
3 P'P 1 
3 P P 1 
BPP 1 


ST OR 4 




tlifSBSssf® 

iffp 


i n r 3 


I B.-1AP ASSEMBLY JB 7 


U f >'> 


Y = 1 

Y S = I R 
X(ISK,74)=CP*(YS-Y) 


I J=( I-* ) *JR+J 

1 FORMAT (5X,8 t 5.7) 

5’* CONTI MU' 

5 0 CO JTI MU - 
R L: TURM 
END 

fortran 

T I 0 N $ U 6 P R C) G R A M REF F: R t. NC \ S 


ES FOR ABOVE ASSEMBLY 


I6LDR — JOB 


PROGRAM UB 7 




